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1 Introduction 

Often signals and system parameters are most conveniently represented as complex-valued vec- 
tors. This occurs, for example, in array processing [1], as well as in communication systems [7] 
when processing narrowband signals using the equivalent complex baseband representation [2]. 
Furthermore, in many important applications one attempts to optimize a scalar rea/-valued mea- 
sure of performance over the complex parameters defining the signal or system of interest. This is 
the case, for example, in LMS adaptive filtering where complex filter coefficients are adapted on 
line. To effect this adaption one attempts to optimize the performance measure by adjustments of 
the coefficients along its stochastic gradient direction [16, 23]. 

However, an often confusing aspect of complex LMS adaptive filtering, and other similar 
gradient-based optimization procedures, is that the partial derivative or gradient used in the adapta- 
tion of complex parameters is not based on the standard complex derivative taught in the standard 
mathematics and engineering complex variables courses [3]-[6], which exists if and only if a func- 
tion of a complex variable z is analytic in z} This is because a nonconstant reaZ- valued function 
of a complex variable is not (complex) analytic and therefore is not differentiable in the standard 
textbook complex-variables sense. 



I.e., comp/ex-analytic. 
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Nonetheless, the same real-valued function viewed as a function of the real-valued real and 
imaginary components of the complex variable can have a (real) gradient when partial derivatives 
are taken with respect to those two (real) components. In this way we can shift from viewing 
the real- valued function as a non-differentiable mapping between C and M to treating it as a dif- 
ferentiable mapping between and R. Indeed, the modern graduate-level textbook in complex 
variables theory by Remmert [12] continually and easily shifts back and forth between the real 
function ^ M or perspective and the complex function C C perspective of a complex or 
real scalar-valued function, 

f{z) = f{r) = f{x,y), 

of a complex variable z = x + j y, 

z eC^r = (^^ eR\ 

In particular, when optimizing a real-valued function of a complex variable z = x + j y one can 
work with the equivalent real gradient of the function viewed as a mapping from to M in lieu 
of a nonexistent complex derivative [14]. However, because the real gradient perspective arises 
within a complex variables framework, a direct reformulation of the problem to the real domain 
is awkward. Instead, it greatly simplifies derivations if one can represent the real gradient as a 
redefined, new complex gradient operator. As we shall see below, the complex gradient is an 
extension of the standard complex derivative to non-complex analytic functions. 

Confusing the issue is the fact that there is no one unique way to consistently define a "complex 
gradient" which applies to (necessarily non-complex-analytic) real-valued functions of a complex 
variable, and authors do not uniformly adhere to the same definition. Thus it is often difficult to 
resolve questions about the nature or derivation of the complex gradient by comparing authors. 
Given the additional fact that typographical errors seem to be rampant these days, it is therefore 
reasonable to be skeptical of the algorithms provided in many textbooks-especially if one is a 
novice in these matters. 

An additional source of confusion arises from the fact that the derivative of a function with 
respect to a vector can be alternatively represented as a row vector or as a column vector when a 
space is Cartesian,^ and both representations can be found in the literature. In this note we carefully 
distinguish between the complex cogradient operator (covariant derivative operator [22]), which 
is a row vector operator, and the associated complex gradient operator which is a vector operator 
which gives the direction of steepest ascent of a real scalar- valued function. 

Because of the constant back-and-forth shift between a real function ("M-calculus") perspective 
and a complex function ("C-calculus") perspective which a careful analysis of nonanalytic complex 
functions requires [12], we refer to the mathematics framework underlying the derivatives given 
in this note as a "CK-calculus." In the following, we start by reviewing some of the properties of 
standard univariate analytic functions, describe the CM-calculus for univariate nonanalytic func- 
tions, and then develop a multivariate second order CK-calculus appropriate for optimizing scalar 
real-valued cost functions of a complex parameter vector. We end the note with some examples. 

^I.e., is Euclidean with identity metric tensor. 
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2 The Derivative of a Holomorphic Function 

Let z = X + jy, for x, y real, denote a complex number and let 

f{z) = u{x,y) + jv{x,y) 

be a general complex- valued function of the complex number z.^ In standard complex variables 
courses it is emphasized that for the complex derivative, 



to exist in a meaningful way it must be independent of the direction with which approaches 
zero in the complex plane. This is a very strong condition to be placed on the function f{z). As 
noted in an introductory comment from the textbook by Flanigan [6] : 

You will leaiTi to appreciate the difference between a complex analytic function (roughly 
a complex-valued function f{z) having a complex derivative f'{z)) and the real functions 
y = f{x) which you differentiated in calculus. Don't be deceived by the similarity of the 
notations f{z), f{x). The complex analytic function f{z) turns out to be much more special, 
enjoying many beautiful properties not shared by the run-of-the-mill function from ordinary 
real calculus. The reason [ • • • ] is that f{x) is merely f{x) whereas the complex analytic 
function f{z) can be written as 



forces a great deal of structure on f{z); moreover, this structure mirrors the structure of the 
harmonic u{x, y) and v{x, y), functions of two real variables.^ 

In particular the following conditions are equivalent statements about a complex function f[z) 
on an open set containing z in the complex plane [6] : 

''Later, in Section 3, we will interchangeably alternate between this notation and the more informative notation 
/(z, z). Other useful representations are f{u, v) and /(a;, y). In this section we look for the (strong) conditions for 
which / : z I— > /(z) € C is differentiable as a mapping C — > C (in which case we say that / is C-differentiable), 
but in subsequent sections we will admit the weaker condition that / : (.t, y) > {u, v) be differentiable as a mapping 
(in which case we say that / is M.-diffe rentiable); see Remmert [12] for a discussion of these different types 
of differentiability. 

'^Quoted from page 2 of reference [6]. Note that in the quote i — V~l whereas in this note we take j — V— T 
following standard electrical engineering practice. 




f{z) = u{x,y) +iv{x,y) 



where z = x + iy and y), f (x, y) are each real- valued harmonic functions related to each 
other in a very strong way: the Cauchy-Riemann equations 



du dv dv du 
dx dy dx dy 



In summary, the deceptively simple hypothesis that 



f'{z) exists 
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• The derivative f'{z) exists and is continuous. 

• The function f{z) is holomorphic (i.e, complex-analytic in z).^ 

• The function f{z) satisfies the Cauchy-Riemann conditions (1). 

• All derivatives of the function f{z) exist and f{z) has a convergent power series. 
Furthermore, it is a simple consequence of the Cauchy-Riemann conditions that 

f{z) = u{x,y) + jv{x,y) 

is holomorphic only if the functions u{x, y) and v{x, y) both satisfy Laplace's equation 

d^u{x,y) d^u{x,y) _ d'^v{x,y) d^v{x,y) _ 

— TT-^ 1 TT^ — — u ana — — — 1 — ^ — — u. 

ox"^ oy^ ox^ oy^ 

Such functions are known as harmonic functions. Thus if either y) or v{x, y) fail to be har- 
monic, the function f{z) is not differentiable.^ 

Although many important complex functions are holomorphic, including the functions 2;", e^, 
ln(2), sin(2), and cos{z), and hence differentiable in the standard complex variables sense, there 
are commonly encountered useful functions which are not: 

• The function f{z) = z, where '2' denotes complex conjugation, fails to satisfy the Cauchy- 
Riemann conditions. 

• The functions f{z) = Re(z) = = x and g{z) = lm(z) = = y fail the Cauchy- 
Riemann conditions. 

• The function f{z) = \zf = zz = x"^ + y'^ is not harmonic. 

^ A function is analytic on some domain if it can be expanded in a convergent power series on that domain. Although 
this condition impHes that the function has derivatives of all orders, analyticity is a stronger condition than infinite 
differentiability as there exist functions which have derivatives of all orders but which cannot be expressed as a power 
series. For a complex-valued function of a complex variable, the term (complex) analytic has been replaced in modern 
mathematics by the entirely synonymous term holomorphic. Thus real-valued power-series-representable functions of 
a real-variable are analytic (real-analytic), while complex-valued power-series-representable functions of a complex- 
variable are holomorphic (complex-analytic). We can now appreciate the merit of distinguishing between holomorphic 
and (real) analytic functions — a function can be nonholomorphic (i.e. non-complex-analytic) in the complex variable 
z = X -\- j y yet still be (real) analytic in the real variables x and y. 

^Because a harmonic function on satisfies the partial differential equation known as Laplace's equation, by 
existence and uniqueness of the solution to this partial differential equation its value is completely determined at 
a point in the interior of any simply connected region which contains that point once the values on the boundary 
(boundary conditions) of that region are specified. This is the reason that contour integration of a complex-analytic 
(holomorphic) function works and that we have the freedom to select that contour to make the integration as easy as 
possible. On the other hand, there is, in general, no equivalent to contour integration for an arbitrary function on M^. 
See the excellent discussion in Flanigan [6]. 
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• Any nonconstant purely real-valued function f{z) (for which it must be the case that v{z, y) = 
0) fails the Cauchy-Riemann condition. In particular the real function f{z) = \z\ = \fYz = 
^/x'^ -\- y'^ is not differentiable^ 

Note in particular, the implication of the above for the problem of minimizing the real- valued 
squared-error loss functional 

i{a) = E{|?7fc - a^fcl^} = E |(?7fc - <fc)(?7fc - =£{6^6^} (2) 

for finite second-order moments stationary scalar complex random variables and 77^., and un- 
known complex constant a = + jay. Using the theory of optimization in Hilbert spaces, the 
minimization can be done by invoking the projection theorem (which is equivalent to the orthogo- 
nality principle) [34]. Alternatively, the minimization can be performed by completing the square. 
Either procedure will result in the Wiener-Hopf equations, which can then be solved for the optimal 
complex coefficient variable a. 

However, if a gradient procedure for determining the optimum is desired, we are immediately 
stymied by the fact that the purely real nonconstant function ^(a) is not complex-analytic (holo- 
morphic) and therefore its derivative with respect to a does not exist in the conventional sense of a 
complex derivative [3]-[6], which applies only to holomorphic functions of a. A way to break this 
impasse will be discussed in the next section. Meanwhile note that all of the real-valued nonholo- 
morphic functions shown above can be viewed as functions of both z and its complex conjugate z, 
as this fact will be of significance in the following discussion. 

3 Extensions of the Complex Derivative - The CM- Calculus 

In this section we continue to focus on functions of a single complex variable z. The primary 
references for the material developed here are Nehari [11], Remmert [12], and Brandwood [14]. 

3.1 A Possible Extension of the Complex Derivative. 

As we have seen, in order for the complex derivative of a function oi z = x + j y, 

f{z) = u{x,y) + jv{x,y), 

to exist in the standard holomorphic sense, the real partial derivatives of u{x, y) and v{x, y) must 
not only exist, they must also satisfy the Cauchy-Riemann conditions (1). As noted by Flanigan 
[6]: "This is much stronger than the mere existence of the partial derivatives." However, the 
"mere existence" of the (real) partial derivatives is necessary and sufficient for a stationary point 

^Thus we have the classic resuh that the only holomorphic real-valued functions are the constant real-valued 
functions. 
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of a (necessarily nonholomorphic) non-constant real-valued functional f{z) to exist when f{z) is 
viewed as a dijferentiable function of the real and imaginary parts of z, i.e., as a function over M^, 

/(z) =/(x,y) :M2^M. (3) 

Thus the trick is to exploit the real vector space structure which underlies C when performing 
gradient-based optimization. In essence, the remainder of this note is concerned with a thorough 
discussion of this "trick." 

Towards this end, it is convenient to define a generalization or extension of the standard partial 
derivative to nonholomorphic functions of z = x + j y that are nonetheless differentiable with 
respect to x and y and which incorporates the real gradient information directly within the complex 
variables framework. After Remmert [12], we will call this the real-derivative, or R-derivative, 
of a possibly nonholomorphic function in order to avoid confusion with the standard complex- 
derivative, or C-derivative, of a holomorphic function which was presented and discussed in the 
previous section. Furthermore, we would like the real-derivative to reduce to the standard complex 
derivative when applied to holomorphic functions. 

Note that if one rewrites the real- valued loss function (2) in terms of purely real quantities, one 
obtains (temporarily suppressing the time dependence, k) 

i[a) = i{a^, ay) = E {e^ + ej} = E {(77^ - a^^^ - ayiyf + {rjy + ay^^ - a^^yf] . (4) 

From this we can easily determine that 



di{ax, a. 



yj 



da^ 

and 



di{ax, a 



yj 



day 



2 E "{Ca^^a; ~l~ ^y^y^ 
2 E {Cx^j/ ^y^x\ 



Together these can be written as 



d , ■ \ „/ \ di{ax,ay) .di{ax,ayj 



da^ day J da^ day 

which looks very similar to the standard result for the real case. 

Indeed, equation (5) is the definition of the generalized complex partial derivative often given in 
engineering textbooks, including references [7]-[9]. However, this is not the definition used in this 
note, which instead follows the formulation presented in [10]-[20]. We do not use the definition 
(5) because it does not reduce to the standard C-derivative for the case when a function /(a) is a 
holomorphic function of the complex variable a. For example, take the simplest case of /(a) = a, 
for which the standard derivative yields ^/(a) = 1- In this case, the definition (5) applied to 
/(a) unfortunately results in the value 0. Thus we will not view the definition (5) as an admissible 
generalization of the standard complex partial derivative, although it does allow the determination 
of the stationary points of i{a).^ 



Hn fact, it is a scaled version of the conjugate M-derivative discussed in the next subsection. 
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3.2 The M-Derivative and Conjugate M-Derivative. 

There are a variety of ways to develop the formalism discussed below (see [11]-[14]). Here, we 
roughly follow the development given in Remmert [12] with additional material drawn from Brand- 
wood [14] and Nehari [11]. 

Note that the nonholomorphic (nonanalytic in the complex variable z) functions given as ex- 
amples in the previous section can all be written in the form /(z, z), where they are holomorphic 
m z = X + jy for fixed z and holomorphic in z = x — jy for fixed z.^ It can be shown that this 
fact is true in general for any complex- or real- valued function 

f{z) = f{z, z) = f{x, y) = u{x, y) + j v{x, y) (6) 

of a complex variable for which the real-valued functions u and v are differentiable as functions 
of the real variables x and y. This fact underlies the development of the so-called Wirtinger cal- 
culus [12] (or, as we shall refer to it later, the CR-calculus.) In essence, the so-called conjugate 
coordinates, 



Conjugate Coordinates: c = (z, z)^ e C x C , z = x + jy and z = x-jy 



(7) 



can serve as a formal substitute for the real r = (x, y)^ representation of the point z = x + jyEC 
[12].^° According to Remmert [12], the calculus of complex variables utilizing this perspective was 
initiated by Henri Poincare (over 100 years ago!) and further developed by Wilhelm Wirtinger in 
the 1920's [10]. Although this methodology has been fruitfully exploited by the German- speaking 
engineering community (see, e.g., references [13] or [31]), it has not generally been appreciated 
by the English speaking engineering community until relatively recently.^' 

For a general complex- or real-valued function /(c) = f{z,z) consider the pair of partial 
derivatives of f (a) formally^^ defined by 



-Derivative of /(c) ^ ^■^^^'^^ 



dz 



df(z z) 

and Conjugate M-Derivative of /(c) = ^ ' ' 



d 

z— const. 



(8) 



^That is, if we make the substitution w = z, they are analytic in w for fixed z, and analytic in z for fixed w. This 
simple insight underlies the development given in Brandwood [14] and Remmert [12]. 

"^Warning! The interchangeable use of the various notational forms of / implicit in the statement f{z) = /(z, z) 
can lead to confusion. To minimize this possibility we define the term "/(z) (z-only)" to mean that f{z) is independent 
of z (and hence is holomorphic) and the term "/(z) {z only)" to mean that f{z) is a function of z only. Otherwise 
there are no restrictions on /(z) ~ f{z, z). 

"An important exception is Brandwood [14] and the work that it has recently influenced such as [1, 15, 16]. 
However, these latter references do not seem to fully appreciate the clarity and ease of computation that the Wirtinger 
calculus (CR-calculus) can provide to the problem of differentiating nonholomorphic function and optimizing real- 
valued functions of complex variables. Perhaps this is do to the fact that [14] did not reference the Wirtinger calculus 
as such, nor cite the rich body of work which had already existed in the mathematics community ([11, 18, 12]). 

'^These statements are formal because one cannot truly vary z — x + j y while keeping z = x — j y constant, and 
vice versa. 
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where the formal partial derivatives are taken to be standard complex partial derivatives (C-derivatives) 
taken with respect to z in the first case and with respect to z in the second. For example, with 

f{z, z) = z^z we have 

— = 2zz and -— = z . 

az oz 

As denoted in (8), we call the first expression the ^-derivative (the real-derivative^ and the second 
expression the conjugate ^.-derivative (or 'R-derivative). 

It is proved in [11, 14, 12] that the R-derivative and iR-derivative formally defined by (8) can 
be equivalently written as^"* 

^ - 1 and (^+ (9) 

dz 2 \dx ^ dy J dz 2 \dx dy J 

where the partial derivatives with respect to x and y are true (i.e., non-formal) partial derivatives of 
the function f(z) = f{x, y), which is always assumed in this note to be differentiable with respect 
to X and y (i.e., to be R-differentiable). Thus it is the right-hand-sides of the expressions given in 
(9) which make rigorous the formal definitions of (8). 

Note that from equation (9) that we immediately have the properties 

?£ = |! = 1 and ^ = ^=0^ (10) 

OZ oz oz oz 

Comments: 

1. The condition || = is true for an M-differentiable function / if and only the Cauchy- 
Riemann conditions are satisfied (see [11, 14, 12]). Thus a function f is holomorphic 
(complex-analytic in z) if and only if it does not depend on the complex conjugated vari- 
able z. I.e., if and only if f{z) = f{z) (z only).^^ 

2. The R-derivative, of an K-differentiable function / is equal to the standard C-derivative, 
f'{z), when f(z, z) is independent of z, i.e., when f{z) = f{z) (z only). 



'^A careful and rigorous analysis of these formal partial derivatives can be found in Remmert [12]. In [12], a 
differentiable complex function / is called C-differentiable while if / is differentiable as a mapping from K^, 
it is said to be real-dijferentiable (R-differentiable) (See Footnote 3). It is shown in [12] that the partial derivatives (8) 
exist if and only if / is R-differentiable. As discussed further below, throughout this note we assume that all functions 
are globally real-analytic (K-analytic), which is a sufficient condition for a function to be globally M-differentiable. 

''^Recall the representation / = /(x, y) = u{x, y) + j v{x, y). Note that the relationships (9) make it clear why the 
partial derivatives (8) exist if and only if / is R-differentiable. (See footnotes 3 and 13). 

'^This obviously provides a simple and powerful characterization of holomorphic and nonholomorphic functions 
and shows the elegance of the Wirtinger calculus formulation based on the use of conjugate coordinates (z, z). Note 
that the two Cauchy-Riemann conditions are replaced by the single condition §j = 0. The reader should reexamine 
the nonholomorphic (nonanalytic in z) functions discussed in the previous section in the light of this condition. 
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3. An M-differentiable function / is holomorphic in z (complex- analytic in z) if and only if 
it does not depend on the variable z, f{z, z) = f{z) (z only), which is true if and only if 

^ = 

To summarize, an E-differentiable function / is holomorphic (complex-analytic in z) if and only 



if f{z) = f{z) (z only), which is true if and only if 



91 

dz 



0, in which case the M-derivative 



coincides with the standard C-derivative, |^ = f'{z). We call the single condition || = the 
Cauchy-Riemann condition for / to be holomorphic: 



Cauchy Riemann Condition: 



dz 



(11) 



Real-Analytic Complex Functions. Throughout the discussion given above we have been mak- 
ing the assumption that a complex function / is real differentiable (M-differentiable). We hence- 
forth make the stronger assumption that complex functions over C are globally real-analytic (R- 
analytic) over M?. As discussed above, and rigorously proved in Remmert [12], R- analytic func- 
tions are M-differentiable and H-differentiable. 

A function f{z) has a power series expansion in the complex variable z, 

f{z) = f{zo) + f{zo){z - zo) + lf"{zo){z - zof + ■■■ + ^/^"n^o)(^ - ZoT + ■■■ 

where the complex coefficient (^o) denotes an n-times C-derivative of f{z) evaluated at the 
point zq, if and only if it is holomorphic in an open neighborhood of zq. If the function f{z) is not 
holomorphic over C, so that the above expansion does not exist, but is nonetheless still M-analytic 
as a mapping from M? to M?, then the real and imaginary parts of f{z) = u{x,y) + jv{x,y), 
z = X + j y, can be expanded in terms of the real variables r = (x, y)^, 

, , , , du{ro) J. d f du{ro)\^ 

u{r) = u{ro) + —^{r - ro) + {r - ro) j (r - Tq) + ■ ■ ■ 

v{r) = v{ro) + —^{r - ro) + [r - ro) \—— j (r - ro) + ■ ■ ■ 

Note that if the R-analytic function is purely real, then f{z) = u{x, y) and we have 

/W = /W + ^(r - r„) + ir - (^)' (r - r„) + • • • 
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Properties of the R- and R-Derivatives. The M-derivative and Tl-derivative are both linear oper- 
ators which obey the product rule of differentiation. The following important and useful properties 
also hold (see references [11, 12]).'^ 

Complex Derivative Identities: 



dz \ dz 



dl ^ fdl 
dz V dz 



(12) 
(13) 



df df 

df = -J-dz + ^dz Differential Rule (14) 



dz dz 
dh(g) dh dg dh dg 

dz dg dz dg dz 



Chain Rule (15) 



dh(g) dh dg dh dg ■ n ■ ..^^ 

= I^IF + TFTF Cham Rule (16) 

dz dg dz dg dz 

As a simple consequence of the above, note that if f{z) is real- valued then / [z] = f{z) so that we 
have the additional very important identity that 



As a simple first application of the above, note that the Tl-derivative of i{a) can be easily 
computed from the definition (2) and the above properties to be 

^^(«) -E{efcefc} = E\^ek + ek^\ = E{0-ek-ekCk} = -Ej^fcCfc} . (18) 



da \ da da ^ 

which is the same result obtained from the "brute force" method based on deriving expanding the 
loss function in terms of the real and imaginary parts of a, followed by computing (5) and then 
using the result (9). Similarly, it can be easily shown that the R-derivative of i{a) is given by 

Note that the results (18) and (19) are the complex conjugates of each other, which is consistent 
with the identity (17). 

We view the pair of formal partial derivatives for a possibly nonholomorphic function defined 
by (8) as the natural generalization of the single complex derivative (C-derivative) of a holomorphic 



'^In the following for z ~ x + j y ws define dz = dx + j dy and dz = dx — j dy, while h{g) = ho g denotes the 
composition of the two function h and g. 
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function. The fact that there are two derivatives under general consideration does not need to be 
developed in elementary standard complex analysis courses where it is usually assumed that / 
is always holomorphic (complex-analytic in z). In the case when / is holomorphic then / is 
independent of z and the conjugate partial derivative is zero, while the extended derivative reduces 
to the standard complex derivative. 

First-Order Optimality Conditions. As mentioned in the introduction, we are often interested 
in optimizing a scalar function with respect to the real and imaginary parts r = {x,y)^ of a 
complex number z = x + jy. It is a standard result from elementary calculus that a first-order 
necessary condition for a point rg = (xq, yo)^ to be an optimum is that this point be a stationary 
point of the loss function. Assuming differentiability, stationarity is equivalent to the condition 
that the partial derivatives of the loss function with respect the parameters r = {x, y)^ vanish at 
the point r = (xq, yo)^- The following fact is an easy consequence of the definitions (8) and is 
discussed in [14]: 

• A necessary and sufficient condition for a real- valued function, f{z) = f{x,y), z = x + j y, 
to have a stationary point with respect to the real parameters r = (x, y)^ G is that its R- 
derivative vanishes. Equivalently, a necessary and sufficient condition for f(z) = f{x, y) 
to have a stationary point with respect to r = (x, yY G is that its M-derivative vanishes. 

For example, setting either of the derivatives ( 1 8) or ( 1 9) to zero results in the so-called Wiener- 
Hopf equations for the optimal MMSE estimate of a. This result can be readily extended to the 
multivariate case, as will be discussed later in this note. 

The Univariate CM-Calculus. As noted in [ 1 2] , the approach we have been describing is known 
as the Wirtinger calculus in the German speaking countries, after the pioneering work of Wilhelm 
Wirtinger in the 1920's [10]. Because this approach is based on being able to apply the calculus 
of real variables to make statements about functions of complex variables, in this note we use the 
term "CR-calculus" interchangeable with "Wirtinger calculus." 

Despite the important insights and ease of computation that it can provide, it is the case that 
the use of conjugate coordinates z and z (which underlies the CE-calculus) is not needed when 
developing the classical univariate theory of holomorphic (complex- analytic in z) functions.'^ It 
is only in the multivariate and/or nonholonomic case that the tools of the CM-calculus begin to 
be indispensable. Therefore it is not developed in the standard courses taught to undergraduate 
engineering and science students in this country [3]-[6] which have changed little in mode of 
presentation from the earliest textbooks. 

''"The differential calculus of these operations ... [is] ... largely irrelevant for classical function theory ..." — 
R. Remmert [12], page 66. 

'^For instance, the widely used textbook by Churchill [3] adheres closely to the format and topics of its first edition 
which was published in 1948. The latest edition (the 7th at the time of this writing) does appear to have one brief 
homework problem on differentiating nonholomorphic functions. 
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Ironically, the elementary textbook by Nehari [11] was an attempt made in 1961 (almost 50 
years ago!) to integrate at least some aspects of the CR-calculus into the elementary treatment of 
functions of a single complex variable.'^ However, because the vast majority of textbooks treat the 
univariate case, as long as the mathematics community, and most of the engineering community, 
was able to avoid dealing with nonholomorphic functions, there was no real need to bring the ideas 
of the CM-calculus into the mainstream univariate textbooks. 

Fortunately, an excellent, sophisticated and extensive introduction to univariate complex vari- 
ables theory and the CK-calculus is available in the textbook by Remmert [12], which is a transla- 
tion from the 1989 German edition. This book also details the historical development of complex 
analysis. The highly recommended Remmert and Nehari texts have been used as primary refer- 
ences for this note (in addition to the papers by Brandwood [14] and, most importantly for the 
second-order analysis given below, van den Bos [25]). 

The Multivariate CM-Calculus. Although one can forgo the tools of the CM-calculus in the 
case of univariate holomorphic functions, this is not the situation in the multivariate holomorphic 
case where mathematicians have long utilized these tools [17]-[20].^° Unfortunately, multivariate 
complex analysis is highly specialized and technically abstruse, and therefore virtually all of the 
standard textbooks are accessible only to the specialist or to the aspiring specialist. It is commonly 
assumed in these textbooks that the reader has great facility with differential geometry, topology, 
calculus on manifolds, and differential forms, in addition to a good grasp of advanced univariate 
complex variables theory. Moreover, because the focus of the theory of multivariate complex 
functions is primarily on holomorphic functions, whereas our concern is the essentially ignored 
(in this literature) case of nonholomorphic real-valued functionals, it appears to be true that only a 
very small part of the material presented in these references is directly useful for our purposes (and 
primarily for creating a rigorous and self-consistent multivariate CK-calculus framework based on 
the results given in the papers by Brandwood [14] and van den Bos [25]). 

The clear presentation by Brandwood [14] provides a highly accessible aspect of the first-order 
multivariate CM-calculus as applied to the problem of finding stationary values for real-valued 
functionals of complex variables. As this is the primary interest of many engineers, this pithy 
paper is a very useful presentation of just those very few theoretical and practical issues which are 
needed to get a clear grasp of the problem. Unfortunately, even twenty years after its publication, 
this paper still is not as widely known as it should be. However, the recent utilization of the Brand- 
wood results in [1, 13, 15, 16] seems to indicate a standardization of the Brandwood presentation 
of the complex gradient into the mainstream textbooks. The results given in the Brandwood pa- 
per [14] are particulary useful when coupled with with the significant extension of Brandwood's 

'''This is still an excellent textbook that is highly recommended for an accessible introduction to the use of deriva- 
tives based on the conjugate coordinates z and z. 

^^"[The CR-calculus] is quite indispensable in the function theory of several variables." — R. Remmert [12], page 
67. 

^'Although, as mentioned in an earlier footnote, Brandwood for some reason did not cite or mention any prior work 
relating to the use of conjugate coordinates or the Wirtinger calculus. 
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results to the problem of computing complex Hessians which has been provided by van den Bos's 
paper [25]. 

At this still relatively early stage in the development of a widely accepted framework for dealing 
with real-valued (nonholomorphic) functions of several complex variables, presumably even the 
increasingly widely used formalism of Brandwood [14] and van den Bos [25] potentially has some 
room for improvement and/or clarification (though this is admittedly a matter of taste). In this 
spirit, and mindful of the increasing acceptance of the approach in [14] and [25], in the remainder 
of this note we develop a multivariate CM-calculus framework that is only slightly different than that 
of [14] and [25], incorporating insights available from the literature on the calculus of multivariate 
complex functions and complex differential manifolds [17]-[20].^^ 

4 Multivariate CM-Calculus 

The remaining sections of this note will provide an expanded discussion and generalized presen- 
tation of the multivariate CR-calculus as presented in Brandwood [14] and van den Bos [25]. The 
discussion given below also utilizes insights gained from references [17, 18, 19, 20, 21, 22]. 

4.1 The Space Z = C". 

We define the n-dimensional column vector z by 

where Zi = Xi + j yi,i = 1, - ■ ■ , n, or, equivalently, 

z = X + j y 

with X = (xi ■ ■ ■ XnY and y = (yi - ■ ■ Vn)^ ■ The space Z = C" is a vector space over the field 
of complex numbers with the standard component-wise definitions of vector addition and scalar 
multiplication. Noting the one-to-one correspondence 

z G C" ^ r = ( ^ ] G 7^ = M^" = X 



it is evident that there exists a natural isomorphism between Z = C" and IZ = R^". 
The conjugate coordinates of z G C" are defined by 

z=(% ■■■ ^„)^gZ = C" 



^^Realistically, one must admit that many, and likely most, practicing engineers will be unlikely to make the move 
from the perspective and tools provided by [14] and [25] (which already enable the engineer to solve most problems 
of practical interest) to that developed in this note, primarily because of the requirement of some familiarity of (or 
willingness to learn) concepts of differential geometry at the level of the earlier chapters of [21] and [22]). 
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We denote the pair of conjugate coordinate vectors (z, z) by 

:^ G c^" = c" X c" 

Noting that c, (z, z), z, (x, y), and r are alternative ways to denote the same point z = x + j y 
in Z = C", for a function 

f : ^ C"" 

throughout this note we will use the convenient (albeit abusive) notation 

f (c) = f (z, z) = f (z) = f (x, y) = f (r) G C"^ 

where z = x + jyGZ = C". We will have more to say about the relationships between these 
representations later on in Section 6 below. 

We further assume that Z = is a Riemannian manifold with a hermitian, positive-definite 
n X n metric tensor fi^ = > 0. This assumption makes every tangent space^^ TzZ = a 
Hilbert space with inner product 

(Vi, V2) = vf fizV2 Vi, V2 G C^. 




4.2 The Cogradient Operator and the Jacobian Matrix 

The Cogradient and Conjugate Cogradient. Define the cogradient and conjugate cogradient 
operators respectively as the row operators^ 



d 

Cogradient Operator: — = 

cz 



dZn I 



(20) 



where (zj, Zi] 
operators are 



Conjugate cogradient Operator: 



dZn / 



(21) 



i = 1, ■ ■ ■ , are conjugate coordinates as discussed earlier and the component 
^-derivatives and R-derivatives defined according to equations (8) and (9), 



d_ 

dzi 



1 
2 



_d_ 

dxi 



J 



d_ 

dyi 



and = - 

dzi 2 



dxi ^ dyi 



(22) 



^^A tangent space at the point z is the space of all differential displacements, dz, at the point z or, alternatively, 
the space of all velocity vectors v = ^ at the point z. These are equivalent statements because dz and v are scaled 
version of each other, dz ~ vdt. The tangent space TzZ = is a linear variety in the space Z = C". Specifically it 
is a copy of C" affinely translated to the point z, = {z} + C". 

-'^The "cogradient" is a covariant operator [22]. It is not itself a gradient, but is the companion to the gradient 
operator defined below. 
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for i 



n. Equivalently, we have 

— -- (—- ■—\ and —-- ( — ■ — 



(23) 



When applying the cogradient operator z is formally treated as a constant, and when 
applying the conjugate cogradient operator z is formally treated as a constant. For example, 
consider the scalar- valued function 

/(c) = /(Z, Z) = Z1Z2 + Z1Z2 . 

For this function we can readily determine by partial differentiation on the Zi and Zi components 
that 



df{c) 



(^2 Zi) and 



df{c) 



{Z2 Zi) 



dz ^ dz 

The Jacobian Matrix. Let f (c) = f (z, z) e C™ be a mapping^^ 

f : Z = C" ^ C"". 

The generalization of the identity (14) yields the vector form of the differential rule,^^ 

Differential Rule 



df{c) = ^dc 



dc 



9f (c) ^ 9f (c) ^_ 
^ -dzA ^dz. 



dz 



dz 



(24) 



where the m x n matrix |^ is called the Jacobian, or Jacobian matrix, of the mapping f , and the 
m X n matrix || the conjugate Jacobian of f . The Jacobian of f is often denoted by Jf and is 
computed by applying the cogradient operator component- wise to f , 



Jf(c) 



(9f(c) 
dz 



oz 



\ dz / 



9. flic) 


9/1 (c) 


dzi 


dzn 


dU(c) 


9/n(c) 


dzi 





(25) 



and similarly the conjugate Jacobian, denoted by Jj is computing by applying the conjugate cogra- 
dient operator component-wise to f , 



gf(c) 
dz 



az 



\ dz / 



/ 8/i(c) 



8/n(c) 



9/1 (c) \ 

dZn * 



a/n(c) 



(26) 



■^^As before the left-hand-sides of (22) and (23) are formal partial derivatives, while the right-hand-sides are actual 
partial derivatives. 

-^It will always be assumed that the components of vector-valued functions are M-differentiable as discussed in 
footnotes (3) and (13). 

^^At this point in our development, the expression ^^^dc only has meaning as a shorthand expression for dz+ 
dz, each term of which must be interpreted formally as z and z cannot be varied independently of each other. 
(Later, we will examine the very special sense in which the a derivative with respect to c itself can make sense.) Also 
note that, unlike the real case discussed in [22], the mapping dz 1-^ di{c) is not linear in dz. Even when interpreted 
formally, the mapping is affine in dz, not linear 
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With this notation we can write the differential rule as 

df (c) = Jf (c) dz + Jf^(c) dz . Differential Rule (27) 
Applying properties (12) and (13) component- wise yields the identities 



9f(c) /9f(c)\ di{c) fdi{c 



Jf (c) and = — ^ = Jf"(c) . (28) 



dz \ dz J dz \ dz 

Note from (28) that, 



However, in the important special case that f (c) is real-valued (in which case f (c) = f(c)) we 
have 

f(c) G ^ Jf(c) = ^ = ^ = Jf^(c). (30) 

dz dz 

With (27) this yields the following important fact which holds for real- valued functions f (c),^^ 



f (c) G ^ rff (c) = Jf (c) dz + Jf(c) rfz = 2Re {Jf (c) dz] . (31) 

Consider the composition of two mappings h : C™ C and g : C" ^ C'", 

h o g = h(g) : C" ^ . 
The vector extensions of the chain rule identities (15) and (16) to h o g are 

9h(g) dh. dg dh dg 



dz dg dz dg dz 

dh{g) dh dg dh dg 

dz dg dz dg dz 

which can be written as 



Chain Rule (32) 
Chain Rule (33) 



Jhog — Jh Jg + Jh (34) 

"^hog ~ Jh Jg + Jh Jg (35) 



^^The real part of a vector (or matrix) is the vector (or matrix) of the real parts. Note that the mapping dz i-^ df{c) 
is not linear 
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Holomorphic Vector- valued Functions. By definition the vector-valued function f(z) is holo- 
morphic (analytic in the complex vector z) if and only if each of its components 



fi{c) = /i(z,z) = fi{zi, 



m 



is holomorphic separately with respect to each of the components Zj,j = 1, ■ ■ ■ ,n. In the refer- 
ences [17, 18, 19, 20] it is shown that f (z) is holomorphic on a domain if and only if it satisfies a 
matrix Cauchy Riemann condition everywhere on the domain: 



di 

Cauchy Riemann Condition: = — = 

oz 



(36) 



This shows that a vector-valued function which is holomorphic on C" must be a function ofz only, 
f (c) = f (z, z) = f (z) (z only). 

Stationary Points of Real- Valued Functionals. Suppose that / is a scalar real-valued function 
from C" to M,29 

f:C-^R; z^/(z). 

As discussed in [14], the first-order differential condition for a real- valued functional / to be 
optimized with respect to the real and imaginary parts of z at the point zq is 



Condition I for a Stationary Point: 



dz 







(37) 



That this fact is true is straightforward to ascertain from equations (20) and (23). An equivalent 
first-order condition for a real-valued functional / to be stationary at the point zq is given by 



Condition II for a Stationary Point: 



dz 







(38) 



The equivalence of the two conditions (37) and (38) is a direct consequence of (28) and the fact 
that / is real-valued. 



Differentiation of Conjugate Coordinates? Note that the use of the notation /(c) as shorthand 
for /(z, z) appears to suggest that it is permissible to take the complex cogradient of /(c) with 
respect to the conjugate coordinates vector c by treating the complex vector c itself as the variable 
of differentiation. This is not correct. Only complex differentiation with respect to the complex 
vectors z and z is well-defined. Thus, from the definition c = col(z, z) G C^", for c viewed as a 
complex 2n-dimensional vector, the correct interpretation of ^/(c) is given by 



^/(z,z) , Jr/(z,z) 
OZ oz 



29' 



The function / is unbolded to indicate its scalar-value status. 
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Thus, for example, we have that 

oc 

which would be true if it were permissible to take the complex cogradient with respect to the 
complex vector c (which it isn't). 

Remarkably, however, below we will show that the 2n-dimensional complex vector c is an 
element of an n-dimensional real vector space and that, as a consequence, it is permissible to take 
the real cogradient with respect to the real vector c! 



Comments. With the machinery developed up to this point, one can solve optimization problems 
which have closed-form solutions to the first-order stationarity conditions. However, to solve 
general nonlinear problems one must often resort to gradient-based iterative methods. Furthermore, 
to verify that the solutions are optimal, one needs to check second order conditions which require 
the construction of the hessian matrix. Therefore, the remainder of this note is primarily concerned 
with the development of the machinery required to construct the gradient and hessian of a scalar- 
valued functional of complex parameters. 



4.3 Biholomorphic Mappings and Change of Coordinates. 

Holomorphic and Biholomorphic Mappings. A vector- valued function f is holomorphic (complex- 
analytic) if its components are holomorphic. In this case the function does not depend on the 
conjugate coordinate z, f (c) = f (z) (z-only), and satisfies the Cauchy-Riemann Condition, 

az 

As a consequence (see (27)), 

9f fz) 

f (z) holomorphic ^ df (z) = Jf (z) rfz = dz . (39) 

Note that when f is holomorphic, the mapping dz i— > di{z) is linear, exactly as in the real case. 
Consider the composition of two mappings h : — > C and g : C" ^ C", 

h o g = h(g) : ^ C^ 

which are both holomorphic. In this case, as a consequence of the Cauchy-Riemann condition 
(36), the second chain rule condition (35) vanishes, J^^g = 0, and the first chain rule condition 
(34) simplifies to 

f and g holomorphic =^ Jhog = JhJg- (40) 
Now consider the holomorphic mapping ^ = f (z). 



d^ = rff(z) = Jf(z)c/z 



(41) 
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and assume that it is invertible, 

z = g{^) = r\^). (42) 

If the invertible function f and its inverse g = f are both holomorphic, then f (equivalently, g) 
is said to be biholomorphic. In this case, we have that 



dT. 



dm 



d$ = U^) d$ = Ji\z) , £ = f(z) 



showing that 



U^) = Ji\z), ^ = f(z). 



(43) 
(44) 



Coordinate Transformations. Admissible coordinates on a space defined over a space of com- 
plex numbers are related via biholomorphic transformations [17, 18, 19, 20]. Thus if z and ^ are 
admissible coordinates on Z = C", there must exist a biholomorphic mapping relating the two 
coordinates, ^ = f (z). This relationship is often denoted in the following (potentially confusing) 



manner. 



dz 



(45) 



z = z(0, dz = ^d^ = U^)d^, ^ = ^^(^) = ^^'(")=(^) (4^^ 

These equations tell us how vectors (elements of any particular tangent space C") properly trans- 
form under a change of coordinates. 

In particular under the change of coordinates ^ = ^(z), a vector v G must transform to its 
new representation w E ^^(^z) according to the 



Vector Transformation Law: w = — v = v 

oz 



(47) 



For the composite coordinate transformation 77 (^(z)), the chain rule yields 



df] dri d$, 

Transformation Chain Rule: tt^ = ttt tt- or Jvo^ = J-n <4 

dz Ok, oz 



(48) 



Finally, applying the chain rule to the cogradient, ||, of a an arbitrary holomorphic function f 



we obtain 



di di dz i if ^ 



This shows that the cogradient, as an operator on holomorphic functions, transforms like 



Cogradient Transformation Law: 



d{-) _d{-) dz _d{-) J 
dz d^ dz ' 



dz ^ 



(49) 
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Note that generally the cogradient transforms quite differently than does a vector. 

Finally the transformation law for the metric tensor under a change of coordinates can be deter- 
mined from the requirement that the inner product must be invariant under a change of coordinates. 
For arbitrary vectors vi, V2 G transformed as 

Wi = Vi G C^(^) 1 = 1,2, 

we have 

(Wi, W2) = wf W2 = vf V2 = vf J^^ V2 = vf V2 = (Vi, V2) . 

This results in the 

(50) 



Metric Tensor Transformation Law: = J^^ J^^ = J, 



5 The Gradient Operator Vz 

l**-Order Approximation of a Real- Valued Function. Let /(c) be a real-valued scalar''^ func- 
tional to be optimized with respect to the real and imaginary parts of the vector z G 2^ = C", 

/ : C" ^ M . 

As a real-valued function, /(c) does not satisfy the Cauchy-Riemann condition (36) and is there- 
fore not holomorphic. 

From (31) we have (with /(z) = /(z, z) = /(c)) that 

ci/(z) = 2Re{J/(z)c/z} = 2Re|^^rfz| . (51) 
This yields the first order relationship 

/(z + dz) = /(z) + 2Re rfzj (52) 

and the corresponding first-order power series approximation 

/(z + Az) ^ /(z) + 2Re 1^ Azj (53) 
which will be rederived by other means in Section 6 below. 



30 



And therefore unbolded. 
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The Complex Gradient of a Real-Valued Function. The relationship (51) defines a nonlinear 
functional, df c{-), on the tangent space C^,^' 



ci/e(v) = 2Re^^v 



V e , c 



z, z 



(54) 



Assuming the existence of a metric tensor ^l^ we can write 



dz 



dz 



H 



H 



(55) 



where Vz/ is the gradient of /, defined as 




(56) 



Consistent with this definition, the gradient operator is defined as 




(57) 



Note the relationships between the gradients and the cogradients. One can show from the coordi- 
nate transformation laws for cogradients and metric tensors that the gradient Vz/ transforms like 
a vector and therefore is a vector, 

Vz/gQ. 



Equations (54) and (55) yield, 

rf/e(v)=2Re{(Vz/, v)} . 

Keeping ||v|| = 1 we want to find the directions v of steepest increase in the value of |(i/c(v)|. We 
have as a consequence of the Cauchy-Schwarz inequality that for all unit vectors v E , 

M/e(v)| = 2 |Re {(Vz/, v)}| < 2 |(Vz/, v)| < 2 ||Vz/|| ||v|| = 2 ||Vz/|| • 

This upper bound is attained if and only if v oc Vz/, showing that the gradient gives the directions 
of steepest increase, with +Vz/ giving the direction of steepest ascent and — Vz/ giving the 
direction of steepest descent. The result (57) is derived in [14] for the special case that the metric 
is Euclidean fi, = /. 



32 



Note that the first-order necessary conditions for a stationary point to exist is given by Vz/ = 0, 

ch does not require kno\ 
/ as is the case in [14]. 



but that it is much easier to apply the simpler condition |f = which does not require knowledge 



of the metric tensor. Of course this distinction vanishes when 



^'Because this operator is nonlinear in dz, unlike the real vector-space case [22], we will avoid calling it a "differ- 
ential operator.". 

^-Therefore one must be careful to ascertain when a result derived in [14] holds in the general case. Also note the 
corresponding notational difference between this note and [14]. We have Vz denoting the gradient operator for the 
general case fl^ I while [14] denotes the gradient operator as for the special case ft^ = I- 
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Comments on Applying the Multivariate CK-Calculus. Because the components of the cogra- 
dient and conjugate cogradient operators (20) and (21) formally behave like partial derivatives 
of functions over real vectors, to use them does not require the development of additional vector 
partial-derivative identities over and above those that already exist for the real vector space case. 
Real vector space identities and procedures for vector partial-differentiation carry over without 
change, provided one first carefully distinguishes between those variables which are to be treated 
like constants and those variables which are to be formally differentiated. 

Thus, although a variety of complex derivative identities are given in various references [14, 15, 
16], there is actually no need to memorize or look up additional "complex derivative identities" 
if one already knows the real derivative identities. In particular, the derivation of the complex 
derivative identities given in references [14, 15, 16] is trivial if one already knows the standard 
real- vector derivative identities. For example, it is obviously the case that 

^(a''z)=a''^ = 0, 



as z is to be treated as a constant when taking partial derivatives with respect to z. Therefore the 
fact that ^a^z = does not have to be memorized as a special complex derivative identity. 

To reiterate, if one already knows the standard gradient identities for real- valued functions of 
real variables, there is no need to memorize additional complex derivative identities?^ Instead, 
one can merely use the regular real derivative identities while keeping track of which complex 
variables are to be treated as constants?'^ This is the approach used to easily derive the complex 
LMS algorithm in the applications section at the end of this note. 

To implement a true gradient descent algorithm, one needs to know the metric tensor. The cor- 
rect gradient, which depends on the metric tensor, is called the "natural gradient" in [24] where it 
is argued that superior performance of gradient descent algorithms in certain statistical parameter 
estimation problems occurs when the natural gradient is used in lieu of the standard "naive" gra- 
dient usually used in such algorithms (where "naive" corresponds to assuming that fi;^ = J even if 
that is not the case). However, the determination of the metric tensor for a specific application can 
be highly nontrivial and the resulting algorithms significantly more complex, as discussed in [24], 
although there are cases where the application of the natural gradient methodology is surprisingly 
straightforward. 

To close this section, we mention that interesting and useful applications of the CM-calculus as 
developed in [14] and [25] can be found in references [13], [26]-[33], and [36], in addition to the 
plentiful material to be found in the textbooks [1], [15], [16], and [23]. 



^^This extra emphasis is made because virtually all of the textbooks (even the exemplary text [15]) provide such 
extended derivative identities and use them to derive results. This sends the message that unless such identities are 
at hand, one cannot solve problems. Also, it places one at the mercy of typographical errors which may occur when 
identities are printed in the textbooks. 

^"^Thus, in the real case, x is the variable to be differentiated in x^^x and we have ^^x^x = 2x"^, while in the 
complex case, if we take z to be treated as constant and z to be the differentiated variable, we have ■^'l^'l = 
z^^z = z^. Note that in both cases we use the differentiation rules for vector differentiation which are developed 
initially for the purely real case once we have decided which variables are to be treated as constant. 
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6 Expansions of a Real-Valued Function on C 



It is common to numerically optimize cost functional using iterative gradient descent-like tech- 
niques. Determination of the gradient of a real- valued loss function via equation (56) allows the 
use of elementary gradient descent optimization, while the linear approximation of a biholomor- 
phic mapping g(^) via (43) enables optimization of the nonlinear least-squares problem using the 
Gauss-Newton algorithm. 

Another commonly used iterative algorithm is the Newton method, which is based on the re- 
peated computation and optimization of the quadratic approximation to the loss function as given 
by a power series expansion to second order. Although the first-order approximation to the loss 
function given by (53) was relatively straight-forward to derive, it is somewhat more work to deter- 
mine the second order approximation, which is the focus of this section and which will be attacked 
using the elegant approach of van den Bos [25].^^ Along the way we will rederive the first order 
approximation (53) and the Hessian matrix of second partial derivatives of a real scalar-valued 
function which is needed to verify the optimality of a solution solving the first order necessary 
conditions. 

6.1 Alternative Coordinate Representations of E = C^. 

Conjugate Coordinate Vectors c G C Form a Real Vector Space. The complex space, C", 
of dimension n naturally has the structure of a real space, M?''\ of dimension 2n, C" ~ M^" , as a 
consequence of the equivalence 



Furthermore, as noted earlier, an alternative representation is given by the set of conjugate 
coordinate vectors 



where C is defined to be the collection of all such vectors c. Note that the set C is obviously a subset 
(and not a vector subspace)^^ of the An dimensional complex vector space C^" . Remarkably, it is 
also a 2n dimensional vector space over the field of real numbers! 

This is straightforward to show. First, in the obvious manner, one can define vector addition 
of any two elements of C. To show closure under scalar multiplication by a real number a is also 
straight forward. 



^^Recall that the Gauss-Newton algorithm is based on iterative re-linearization of a nonlinear model z w g(^)- 
^^A detailed exposition of the second order case is given by Abatzoglou, Mendel, & Harada in [36]. See also 

[32]. The references [36], [25] and [32] all develop the complex Newton algorithm, although with somewhat different 

notation. 

^^It is, in fact, a 27i dimensional submanifold of the space C^" w R*". 




2n 
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Note that this homogeneity property obviously fails when a is complex. 

To demonstrate that C is 2n dimensional, we will construct below the one-to-one transforma- 
tion, J, which maps C onto TZ, and vice versa, thereby showing that C and TZ are isomorphic, 
C ~ 7^. In this manner C and TZ are shown to be alternative, but entirely equivalent (including 
their dimensions), real coordinate representations for Z = V. The coordinate transformation J is 
a linear mapping, and therefore also corresponds to the Jacobian of the transformation between the 
coordinate system IZ and the coordinate system C. 

In summary, we have available three vector space coordinate representations for representing 
complex vectors z = x + j y. The first is the canonical n-dimensional vector space of complex 
vectors z G 2 = C" itself. The second is the canonical 2n-dimensional real vector space of vectors 
r = col(x, y) ElZ = which arises from the natural correspondence C" ~ M^". The third is 
the 2n-dimensional real vector space of vectors c G C C C^", C ~ M^". 

Because C can be alternatively viewed as a complex subset of C^" or as a real vector space iso- 
morphic to M^", we actually have a fourth "representation"; namely the non-vector space complex- 
vector perspective of elements of C as elements of the space C^", c = col(z, z).^^ This perspective 
is just the (z, z) perspective used above to analyze general, possibly nonholomorphic, functions 
/(z) = /(z,z). 

In order to avoid confusion, we will refer to these two alternative interpretations of c G C C 
C^" as the c-real case (respectively, the C-real case) for when we consider the vector c G C ~ M^" 
(respectively, the real vector space C ^ R^"), and the c-complex case (respectively, the C-complex 
case) when we consider a vector c G C C C^" (respectively, the complex subset C C C^").^^ These 
two different perspectives of C are used throughout the remainder of this note. 

Coordinate Transformations and Jacobians. From the fact that 

z = X + j y and z = x — j y 

it is easily shown that 

where / is the n x n identity matrix. Defining"^" 

^^Since when viewed as a subset of C^" the set C is not a subspace, this view of C does not resuh in a true coordinate 
representation. 

^^In the latter case c = col(z, z) is understood in terms of the behavior and properties of its components, especially 
for differentiation purposes because, as mentioned earlier, in the complex case the derivative ^ is not well-defined in 
itself, but is defined in terms of the formal derivatives with respect to z and z. As we shall discover below, in the c-real 
case, the derivative ^ is a true real derivative which is well understood in terms of the behavior of the derivative 

'^''Except for a trivial reordering of the elements of r = (x^ y^)^, this is the transformation proposed and utilized 
by van den Bos [25], who claims in [29] to have been inspired to do so by Remmert. (See, e.g., the discussion on page 
87 of [12].) 
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(58) 



then results in the mapping 
It is easily determined that 

so that we have the inverse mapping 



c(r) = J r . (59) 



J-^ = (60) 



r = r(c) = J-^c = ^J^c. (61) 

Because the mapping between IZ and C is linear, one-to-one, and onto, both of these spaces are 
obviously isomorphic real vector spaces of dimension 2n. The mappings (59) and (61) therefore 
correspond to an admissible coordinate transformation between the c and r representations of 
z G 2. Consistent with this fact, we henceforth assume that the real vector calculus (including all 
of the vector derivative identities) apply to real-valued functions over C. 

Note that for the coordinate transformation c = c(r) = Jr we have the Jacobian 

Je^^c(r) = |-Jr = J (62) 
or or 

showing that J is also the Jacobian of the coordinate transformation from TZ to C.^^ The Jacobian 
of the inverse transformation r = r(c) is given by 

Jr = = = V. (63) 

Of course, then, we have the differential relationships 

dc dv 1 

dc = — dr = Jc dr = Jdr and dr = — dc = Jj-dc = -J^dc (64) 
dr dc 2 



which correspond to the first-order relationships 



42 



1 st-Order Relationships: Ac = Jc Ar = J Ar and Ar = Jr Ac = ^J^Ac 



(65) 



where the Jacobian J is given by (60) and 



Ac=(^^) and Ar=(f;) (66) 



'^'We have just proved, of course, the general property of Unear operators that they are their own Jacobians. 

'^^For a general, nonlinear, coordinate transformation these finite-difference (non-infinitesimal) first-order relation- 
ships would be approximate. However, because the coordinate transformation considered here happens to be linear, 
the relationships are exact. 
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The Cogradient with respect to the Real Conjugate Coordinates Vector c. The reader might 
well wonder why we didn't just point out that (64) and (65) are merely simple consequences of 
the linear nature of the coordinate transformations (59) and (61), and thereby skip the interme- 
diate steps given above. The point is that once we have identified the Jacobian of a coordinate 
transformation over a real manifold, we can readily transform between different coordinate rep- 
resentations of all vector-like (contravariant) objects, such as the gradient of a functional, and 
between all covector-like (covariant) objects, such as the cogradient of a functional, over that 
manifold. Indeed, as a consequence of this fact we immediately have the important cogradient 
operator transformations 



Cogradient Transf s: = ^ j, = r and ^ = ^ = ^ j 

oc or 2 or or Oc oc 



(67) 



with the Jacobian J given by (58) and Jr = J^.^- 

Equation (67) is very important as it allows us to easily, yet rigorously, define the cogradient 
taken with respect to c as a true (nonformal) differential operator provided that we view c as an 

m 

dc 

QnH tVxi=^ "r\nllV\Qr»lr" trcinc^rwmcitirwi 

dr 



element of the real coordinate representation space C. The cogradient ^ is well-defined in terms 
of the cogradient ^ and the "puUback" transformation. 



dc 2 dr 

dc 



This shows that which was originally defined in terms of the cogradient and conjugate cogra- 
dients taken with respect to z (the c-complex interpretation of ^), can be treated as a real differ- 
ential operator with respect to the "real" vector c (the c-real interpretation of ^)-^^ 

Complex Conjugation. It is easily determined that the operation of complex conjugation, z 
z, is a nonlinear mapping on Z = C" . Consider a general element C G C^" written as 

Ctop \ ^^2n_^n^ ^-^^^ ^^^^ ^ ^^^^^^^ 



C bottom 

2n 



Of course the operation of complex conjugation on C \ C, is, in general, a nonlinear mapping. 

Now consider the linear operation of swapping the top and bottom elements of C C — C, 
defined as 

^ _ I Ctop j ^ _ I Cbottom j — ( ^1 ( ^'°P ) ~ ^C 

Cbottom/ V Ctop / VCbottom/ 



'^^Thus we can directly differentiate an expression like c^ilc with respect to c using the standard identities of real 
vector calculus. (The fact that these identities hold for the r calculus and be used to prove their validity for the c-real 
calculus.) More problematic is an expression like c^ilc. It is not appropriate to take the complex derivative of this 
expression with respect to the complex vector c because c, as an element of C" is subject to constraints amongst its 
components. Instead (see immediately below) one can use the identity c = c = 5*0 to obtain c^ilc = Silc which 
can then be differentiated with respect to c. Of course, this latter approach can fail if Silc cannot be interpreted 
in some appropriate sense in the field of real numbers. Note that real versus complex differentiation of c^Qc with 
respect to c would differ by a factor of 2. 
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where 

is the swap operator on C^" which obeys the properties 

s = s^ = s-^, 

showing that S is symmetric and its own inverse, S"^ = I. Note that, in general, swapping is not 
equal to complex conjugation, C 7^ C 

The swap operator S will be used extensively throughout the remainder of this note, so it is 
important to become comfortable with its use and manipulation. The swap operator is a block 
permutation matrix which permutes (swaps)"^"^ blocks of rows or blocks of columns depending on 
whether 5* premultiplies or postmultiplies a matrix. Specifically, let a 2n x 2?7, matrix A be block 
partitioned as 

'All A12 

A2I A22^ 

Then premultiplication by S results in a block swap of the top n rows en masse with the bottom n 

rows,'*^ 

A21 A22 
All Ai2^ 

Alternatively, postmultiplication by S results in a block swap of the first n columns with the last n 
columns 

%2 ^11' 
A22 A2l^ 

It is also useful to note the result of a "sandwiching" by S, 

A22 A21 

A12 An. 



A 



SA 



AS 



SAS = A 



Because S permutes n rows (or columns), it is a product of n elementary permutation matrices, 
each of which is known to have a determinant which evaluates to —1. As an easy consequence of 
this, we have 

det5 = (-1)". 

Other important properties of the swap operator S will be developed as we proceed. 

Now note that the subset C E C^" contains precisely those elements of C^" for which the 
operations of swapping and complex conjugation coincide. 



c = {c G c^" C = c} c c 



^2n 



^'^"Permutation" is just a fancy term for "swapping." 

^^Matrix premultiplication of A by any matrix always yields a row operation. 

46a 



Matrix postmultiplication of A by any matrix always yields a column operation. The fact that pre- and post- 
multiplication yield different actions on A is an interesting and illuminating way to interpret the fact that matrix 
multiplication is noncommutative, MA ^ AM. 
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and thus it is true by construction that c E C obeys c = c, even though swapping and complex 
conjugation are different operations on C^". Now although C is not a subspace of the complex 
vector space C^", it is a real vector space in its own right. We see that the linear operation of 
component swapping on the C-space coordinate representation of Z = C" is exactly equivalent 
to the nonlinear operation of complex conjugation on Z. It is important to note that complex 
conjugation and coordinate swapping represent different operations on a vector c when c is viewed 
as an element of C^"/^ 

We can view the linear swap mapping S : C ^ C as a coordinate transformation (a coordinate 
"reparameterization"), c = c = Sc, on C. Because S is linear, the Jacobian of this transformation 
is just S itself. Thus from the cogradient transformation property we obtain the useful identity 

oc oc oc 

It is also straightforward to show that 

/ = h^S^ (69) 

for J given by (58) 

Let us now turn to the alternative coordinate representation given by vectors r in the space TZ = 
M^". Specifically, consider the TZ coordinate vector r corresponding to the change of coordinates 
r = I J^c. Since the vector r is real, it is its own complex conjugate, r = r."^*^ Complex conjugation 
of z is the nonlinear mapping in C"^ 

z = x + jy^z = x + j (-y) , 

and corresponds in the representation space TZ to the linear mapping^^ 

:)-K->(o-°)e)-^ 



where C is the conjugation matrix 



Note that 



C={1 ^j]. (70) 



c = c^ = 



1 

"^^As mentioned earlier, c, in a sense, does "double duty" as a representation for z; once as a (true coordinate) 
representation of z in the real vector space C, and alternatively as a "representation" of z in the "doubled up" complex 
space C'^" = C" x C". In the development given below, we will switch between these two perspectives of c. 

"^^Note that our theoretical developments are consistent with this requirement, as 

r = -JWe) - ij^c = ij^c = ij^5c = ij^5Jr = /r = r . 
2^ ' 2 2 2 2 

■^^We refer to f as "r-check." 
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i.e., that C is symmetric, C = C , and its own inverse, C = I. It is straightforward to show that 

C = ij^5J (71) 
2 

which can be compared to (69). Finally, it is straightforward to show that 

c = Jr^c = c = Jf. (72) 

To summarize, we can represent the complex vector z by either c or r, where c has two inter- 
pretations (as a complex vector, "c-complex", in C^", or as an element, "c-real", of the real vector 
space C ~ M^"), and we can represent the complex conjugate z by c, c, or f . And complex conju- 
gation, which is a nonlinear operation in C", corresponds to linear operators in the 2n-dimensional 
isomorphic real vector spaces C and TZ. 

6.2 Low Order Series Expansions of a Real-Valued Scalar Function. 

By noting that a real-valued scalar function of complex variables can be viewed as a function of 
either r or c-real or c-complex or z, 

/(r) = /(c) = /(z), 

it is evident that one should be able to represent / as a power series in any of these representations. 
Following the line of attack pursued by van den Bos in [25], by exploiting the relationships (65) 
and (67) we will readily show the equivalence up to second order in a power series expansion of /. 

Up to second order, the multivariate power series expansion of the real-valued function / 
viewed as an analytic function of vector r G 7^ is given as. 



df(r) 1 

2nd-0rder Expansion in r: /(r + Ar) = /(r) + \- ' Ar + - Ar^ Hrr(r) Ar + h.o.t. 

or 2 



where^° 



(73) 



^rr(p) = ^ (^-^J for p,re7^ (74) 

is the real r-Hessian matrix of second partial derivatives of the real- valued function /(r) with 
respect to the components of r. It is well known that a real Hessian is symmetric, 

/T 



rr 



However, there is no general guarantee that the Hessian will be a positive definite or positive 
semidefinite matrix. 

It is assumed that the terms /(r) and /(r + Ar) be readily expressed in terms of c and c + Ac 
or z and z + Az. Our goal is to determine the proper expression of the linear and quadratic terms 
of (73) in terms of c and Ac or z and Az. 



^°When no confusion can arise, one usually drops the subscripts on the Hessian and uses the simpler notation 
7i(p) = Ti-rrip)- Note that the Hessian is the matrix of second partial derivatives of a real-valued scalar function. 
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Scalar Products and Quadratic Forms on the Real Vector Space C. Consider two vectors 
c = col(z, z) e C and s = col(^,^) G C. The scalar product for any two such vectors in C-real 
(i.e., in the real vector space C ~ M^") is defined by 

(c, s) ^c^Ss = cFs = c^s = z^^ + z^l = z^l + = 2 Re z^| . 

The row vector c^S = is a linear functional which maps the elements of C-real into the real 
numbers. The set of all such linear functionals is a vector space itself and is known as the dual 
space, C*,ofC [34, 35]. The elements of C* are known as dual vectors or covectors, and the terms 
"dual vector", "co vector", and "linear functional" should all be taken to be synonymous. Given a 
vector c G C, there is a natural one-to-one mapping between c and a corresponding dual vector, c* 
in C* defined by^^ 

Henceforth it is understood that scalar-product expressions like 

a'^s or c^h 

where s G C and c G C are known to be elements of C are only meaningful if a and b are also 
elements ofC. Thus, it must be the case that both vectors in a scalar product must belong to C if it 
is the case that one of them does, otherwise we view the resulting numerical value as nonsensical. 

Thus, for a real-valued function of up to quadratic order in a vector c G C, 

/(c) = a + b^c + ^c^Afc = a + b^c+ ^c^s, s = Mc, (75) 

to be well-posed, it must be the case that a G M, b G C,^^ and s = Mc G CP Thus, as we proceed 
to derive various first and second order functions of the form (75), we will need to check for these 
conditions. If the conditions are met, we will say that vector b and the operator M; the terms b^c 
and c^Mc; and the entire quadratic form itself, are admissible (or meaningful). 

Thus b is admissible if and only if b G C, and M is admissible if and only if M is a linear 
mapping from C to C, M G £(C, C). 

To test whether a vector b G C^" belongs to C is straightforward: 

heC ^h = Sh. (76) 

It is somewhat more work to develop a test to determine if a matrix M G C^"^^" has the 
property that it is a linear mapping from C to C, 

M G £(C, C) = {M I Mc G C, Vc G C and M is linear } c £(C^", C^") = C^"^^". 

^'Warning! Do not confuse the dual vector (linear functional) c* with an adjoint operator, which is often also 
denoted using the "star" notation. 

^-I.e., that be a bona fide linear functional on C, = b* e C* . 

^^I.e., because = c* G C* , is a linear functional on C, it must have a legitimate object s to operate on, namely 
an element s = Mc <E C. 
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Note that the fact that £(C, C) C £(C^", C^") is just the statement that any matrix which maps 
from C C C^" to C C C^" is also obviously a linear mapping from C^" to C^". However, this is 
just a subset statement; it is not a subspace statement. This is because £(C, C) is a real vector space 
of linear operators,^"^ while £(C^", C^") is a complex vector space of linear operators. Because 
they are vector spaces over different fields, they cannot have a vector-subspace/vector-parent-space 
relationship to each other. 

To determine necessary and sufficient conditions for a matrix M G C^"^^" to be an element 
of C{C, C) suppose that the vector c = col(z, z) E C always maps to a vector s = col(^, ^) E C 
under the action of M, s = Mc. Expressed in block matrix form, this relationship is 

fC] ^ (Mn Mu\ (A 
\M2i M22J \zj- 

The first block row of this matrix equation yields the conditions 

I = Mnz + M12Z 
while the complex conjugate of the second block row yields 

$, = M22Z + M21Z 

and subtracting these two sets of equations results in the following condition on the block elements 
of M, 

(Mn - M22)z + (M12 - M2i)z = . 
With z = X + j y, this splits into the two sets of conditions, 

[(Mn - M22) + (M12 - M2i)]x = 

and 

[(Mn - M22) - (M12 - M2i)]y = 0. 
Since these equations must hold for any x and y, they are equivalent to 

(Mn - M22) + (M12 - M21) = 

and 

(Mn - M22) - (M12 - M21) = 0. 

Finally, adding and subtracting these two equations yields the necessary and sufficient conditions 
for M to admissible (i.e., to be a mapping from C to C), 

M = ( I^^M E C^"^^" is an element of /:(C,C) iff Mn = M22 and M12 = M21 • (77) 

\M21 M22 / 

^"^I.e., a vector space over the field of real numbers. 
^^I.e., a vector space over the field of complex numbers. 
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This necessary and sufficient admissibility condition is more conveniently expressed in the follow- 
ing equivalent form, 

M e CiC, C)^M = SMS ^M = SMS (78) 
which is straightforward to verify. 

Given an arbitrary matrix M E C^"^^", we can define a natural mapping of M into C(C,C) C 

Q2nx2n 

P(M)4^1±^e£(C.C). (79) 
in which case the admissibility condition (78) has an equivalent restatement as 

M e £(C,C) ^ P(M) = M. (80) 
It is straightforward to demonstrate that 

VM G C^'^x^", P(P(M)) = P(M) . (81) 

I.e., P is an idempotent mapping of C^"^^" onto C{C, C), P^ = P. However, as things currently 
stand P is not a linear operator (the action of complex conjugation precludes this) nor a projection 
operator in the conventional sense of projecting onto a lower dimensional subspace as its range 
space is not a subspace of its domain space. (However, with some additional work, one can rea- 
sonably interpret P as a projector of the manifold C^" onto the submanifold C C C^" in some 
sense.^^) 

A final important fact is that if M e C^"''^'' is invertible, then M G C{C,C) if and only if 
E C{C, C), which we state equivalently as 

Let M be invertible, then P(M) = M iff P(M-^) = M'K (82) 

I.e., if an invertible matrix M is admissible, then M~^ is admissible. The proof is straightforward: 

M = SMS and M invertible 

^M~^ = {SMSy^ 

= s{My^s 

= SJF^S. 



s'-With C2"x2,i ^ ]g4«x4n ^ j^i6n= and C{C,C)^ /;(R2",M2"-) ^ ]^2nx2,i ^ ]^4n=^ j( reasonable to view P as 
a linear projection operator from the real vector space M^^"" onto the real vector subspace R^" of R^". This allows 
us to interpret P as a projection operator from the manifold C^" onto the submanifold C C C^". Once we know that 
P is a linear mapping from C^" into C^", we can then compute its adjoint operator, P*, and then test to see if its 
self-adjoint. If it is, then the projection operator P is, in fact, an orthogonal projection operator. 
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First Order Expansions. Up to first order, tlie power series expansion of the real- valued function 
/ viewed as a function of r G 7^ is 



df(r) 

First-Order Expansion in r: /(r + Ar) = /(r) H — Ar + h.o.t. 

OT 



(83) 



Focussing our attention first on the linear term Ar, and using the c-real vector space 
interpretation of c, namely that c G C where, as discussed above, C is a 2?i-dimensional coordinate 
space isomorphic to M^", we have 

df df 

— Ar = — J^^ Ac (from equation (65)) 
or or 

df 

= —- Ac (from equation (67)) 

oc 

which yields the first order expansion of / in terms of the parameterization in c. 



df(c) 

First-Order Expansion in c: /(c + Ac) = /(c) H — Ac + h.o.t. 

oc 



(84) 



Note that Ac is real valued. Furthermore, as a consequence of the fact that with /(c) real- 



valued we have 



^ dc J \ dc J \ dc 

the quantity ( j satisfies the necessary and sufficient condition given in (76) that 

dc 

Thus G C* and the term Ac is admissible in the sense defined earlier. Note that an 

0c ac 

equivalent condition for the term Ac to be admissible is that 

ac 



which is true if and only if 



dc 



G C. 



This shows a simple inspection of itself can be performed to test for admissibility of the 
first-order term.^^ 



57 



In this note, the first order expansion (84) is doing double duty in that it is simultaneously standing for the c-real 
expansion and the c-complex expansion. A more careful development would make this distinction explicit, in which 

case one would more carefully explore the distinction between ( ) versus ( ) in the first-order term. 



Because this note has already become rather notationally tedious, this option for greater precision has been declined. 
However, greater care must therefore be made when switching between the C-real and C-complex perspectives. 
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As discussed above, to be meaningful as a true derivative, the derivative with respect to c has to 
be interpreted as a real derivative. This is provided by the c-real interpretation of (84). In addition, 
(84) has a c-complex interpretation for which the partial derivative with respect to c is not well- 
defined as a complex derivative as it stands, but rather only makes sense as a shorthand notation 
for simultaneously taking the complex derivatives with respect to z and z, 

d__ fd_ d_ 
dc \dz' dz 

Thus, to work in the domain of complex derivatives, we must move to the c-complex perspective 
c = col(z, z), and then break c apart so that we can work with expressions explicitly involving z 
and z, exploiting the fact that the formal partial derivatives with respect to z and z are well defined. 

Noting that 

i) and Ac^(^^) 



we obtain 



-^Ac = -LAz + -^Az 
oc oz oz 

= TT- Az + Az 
Oz Oz 

= 2Re|^Az 
[ oz 



(/is real- valued) 



which yields the first-order expansion of / in terms of the parameterization in z. 



First-Order Expansion in z: /(z + Az) = /(z) + 2 Re 



9f A 

oz 



h.o.t. 



(85) 



This is the rederivation of (53) promised earlier. Note that (85) makes explicit the relationship 
which is implied in the c-complex interpretation of (84). 

We also summarize our intermediate results concerning the linear term in a power series ex- 
pansion using the r, c or z representations. 



Linear- Term Relationships: 



dr 



Ar 



dc 



Ac = 2 Re 



dz 



Az 



(86) 



The derivative in the first expression is a real derivative. The derivative in the second expression 
is interpreted as a real derivative (the c-real interpretation). The derivative in the last expression 
is a complex derivative; it corresponds to the c-complex interpretation of the second term in (86). 
Note that all of the linear terms are real valued. 

We now have determined the first-order expansion of / in terms of r, c, and z. To construct 
the second-order expansion it remains to examine the second-order term in (73) and some of the 
properties of the real Hessian matrix (74) which completely specifies that term. 
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Second Order Expansions. Note from (73) that knowledge of the real Hessian matrix Hrr com- 
pletely specifies the second order term in the real power series expansion of / with respect to r. 
The goal which naturally presents itself to us at this point is now to reexpress this quadratic-order 
term in terms of c, which we indeed proceed to do. However, because the canonical coordinates 
vector c has two interpretations, one as a shorthand for the pair (z, z) (the c-complex perspective) 
and the other as an element of a real vector space (the c-real perspective), we will rewrite the sec- 
ond order term in two different forms, one (the c-complex form) involving the c-complex Hessian 
matrix 

^ d /dfjvy 

dc V dc 



KM 



H 



for v,ceC C C^' 



2n 



(87) 



and the other (the c-real form) involving the c-real Hessian matrix 



d fdfiv) 



dc V dc 



T 



for v,c E C 



D2n 



(88) 



In (87), the derivative with respect to c only has meaning as a short-hand for , ^) . In (88), the 
derivative with respect to c is well-defined via the c-real interpretation. 

It is straightforward to show a relationship between the real Hessian TYrr and the c-complex 
Hessian TY^^., 




(from equation (67)) 



(from equation (67)) 



The resulting important relationship 



C 

CC 



(89) 



between the real and c-complex Hessians was derived in [25] based on the there unjustified (but 
true) assumption that the second-order terms of the powers series expansions of / in terms of r 
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and c-complex must be equal. Here, we reverse this order of reasoning, and will show below the 
equality of the second order terms in the c-complex and r expansions as a consequence of (89). 

Note from (60) that 

^cc = ^J^rrJ''. (90) 

Recalling that the Hessian Hrr is a symmetric matrix,^^ it is evident from (90) that Ti^^, is Hermi- 

tian^^ 

= (Kef 

(and hence, like Hrr, has real eigenvalues), and positive definite (semidefinite) if and only Hrr is 
positive definite (semidefinite). 

As noted by van den Bos [25], one can now readily relate the values of the eigenvalues of Ti^^ 
and Hrr from the fact, which follows from (60) and (90), that 

nl-XI = \j Hrr - -JJ^ = -J {Hrr " 2X1) . 

4 2 4 

This shows that the eigenvalues of the real Hessian matrix are twice the size of the eigenvalues of 
the complex Hessian matrix (and, as a consequence, must share the same condition number).^" 

Focussing our attention now on the second order term of (73), we have 
^Ar^T^rrAr = ^Ar^T^rrAr 

= ^Ar^J^H^JAr (From equation (89)) 
= ^ Ac^ Ac , (From equation (65)) 

thereby showing the equality of the second order terms in an expansion of a real- valued function / 
either in terms of r or c-complex,^^ 

^Ar^T^rr Ar = ^Ac^K'c Ac. (91) 

Note that both of these terms are real valued. 

With the proof of the equalities 86 and 91, we have (almost) completed a derivation of the 



df(c) 1 

2nd-Order Expansion in c-Complex: /(c + Ac) = /(c) + i; ' Ac + -Ac^ Hcc(c) Ac + h.o.t. 

oc 2 



(92) 



^^In the real case, this is a general property of the matrix of second partial derivatives of a scalar function. 

H 

59/ 



^^As expected, as this is a general property of the matrix of partial derivatives (^"g^) of ^i^Y real-valued 
function /(z). 

^"For a Hermitian matrix, the singular values are the absolute values of the (real) eigenvalues. Therefore the condi- 
tion number, which is the ratio of the largest to the smallest eigenvalue (assuming a full rank matrix) is given by the 
ratio of the largest to smallest eigenvalue magnitude. 

*'And thereby providing a proof of this assumed equality in [25]. 
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where the c-complex Hessian H^^. is given by equation (87) and is related to the real hessian Hrr 
by equations (89) and (90). Note that all of the terms in (92) are real valued. The derivation has not 
been fully completed because we have not verified that Ac^ ^cc(c) admissible in the sense 
defined above. The derivation will be fully completed once we have verified that TY^^, G £(C, C), 
which we will do below. 

The c-complex expansion (92) is not differentiable with respect to c-complex itself, which is 
not well defined, but, if differentiation is required, should be instead interpreted as a short-hand, 
or implicit, statement involving z and z, for which derivatives are well defined. To explicitly show 
the second order expansion of the real-valued function / in terms of the complex vectors z and z, 
it is convenient to define the quantities 

Ti-zz = TT ( IT ] ' ^^^ = 7F iir ] ' = 15^ > and Hgz = 7p • (93) 

oz \oz J oz \oz J oz \oz J oz \oz J 

With ^ = , we also have from (87) and the definitions (93) that 

y ''■zz ''■zz J 

Thus, using the earlier proven property that TY^c is Hermitian, TL^^ = {H^^)^ , we immediately 
have from (94) the Hermitian conjugate conditions 

n.z = ng and n^z = ng (95) 

which also hold for z and z replaced by z and z respectively. 

Some additional useful properties can be shown to be true for the block components of (94) de- 
fined in (93). First note that as a consequence of / being a real- valued function, it is straightforward 
to show the validity of the conjugation conditions 



' '-CC ' '-cc 



or, equivalently. 



^zz = ^zz and TYzz = TY^z , (96) 
which also hold for z and z replaced by z and z respectively. It is also straightforward to show that 



^CC ~ S^CC^ — ^ ^CC ^ 5 



for S = S'^ = (showing that H^^ and are related by a similarity transformation and there- 
fore share the same eigenvalues^^), which is precisely the necessary and sufficient condition (78) 
that the Hessian matrix Ti^^ is admissible, H^c ^ ^{C, C). This verifies that the term Ac^H'^^Ac 



^^Their eigenvectors are complex conjugates of each other, as reflected in the similarity transformation being given 
by the swap operator S 
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is admissible and provides the completion of the proof of the validity of (92) promised earlier. 
Finally, note that properties (96) and (95) yield the conjugate symmetry conditions, 

71.. = Hl^ and H.^ = Hl^, (97) 

which also hold for z and z replaced by z and z respectively. 

From equations (66), (91), and (94) we can now expand the second order term in (73) as follows 

-Ar^T^rrAr = -Ac" Hi, Ac 

= i (Az^KzAz + Az^T^^^Az + Az^KzAz + Az^T^^^Az) 
= Re { Az^Hzz Az + Az^Ti^z Az } 

where the last step follows as a consequence of (96).^^ Thus, we have so-far determined that 

^Ar^ rtrr Ar = ^Ac^ Hi, Ac = Re { Az^H^zAz + Az^H^zAz} . (98) 

Combining the results given in (73), (86), and (98) yields the desired expression for the second 
order expansion of / in terms of z. 



2"''-OrderExp.inz: /(z + Az) = /(z) + 2 Re 1 1^ Azj + Re { Az^T^^^Az + Az^HzzAz} + h.o.t. 



(99) 

We note in passing that Equation (99) is exactly the same expression given as Equation (A.7) 
of reference [36] and Equation (8) of reference [32], which were both derived via an alternative 
procedure. 

The c-complex expansion shown in Equation (92) is one of two possible alternative second- 
order representations in c for /(c) (the other being the c-real expansion), and was used as the 
starting point of the theoretical developments leading to the z-expansion (99). We now turn to the 
development of the c-real expansion of /(c), which will be accomplished by writing the second 
order term of the quadratic expansion in terms of the c-real Hessian H^,. 

From the definitions (88), (87), and (93), and using the fact that ^ = , it is straight- 
forward to show that 

njR [Ti-zz ^zz \ c I T^zz 7i,zz\ /'inA\ 

^cc = \n-i = I -ly _ -ly,, (lOU) 



or^^ 



^cc = ^cc = 'S'TYcc = Til-^S. (101) 



^^Alternatively, the last step also follows as a consequence of (95) 
Alternative derivations are possible. For example, H^^. = = ^(Ie)^ = ^[^^ 

(i)^ = ^1; (i)"^ ^^cc ^ W?c = Snl^^ noting that S = = S'^ 



dc \ dc I ~ dc [ dc ~ ac \ ac ^ 
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Note from the first equality in (100) and the conjugate symmetry conditions (97) that the c-real 
Hessian is symmetric 

K^ = iKy. (102) 

Equivalently, 

sni, = {sniy. (103) 

Let the Singular Value Decomposition (SVD) of TY^c be 

then from (101) the SVD oiUl^ is given by 

= U'EV^ , U' = SU 

showing that Ti^,, and Ti^^. share the same singular values, and hence the same condition number 
(which is given by the ratio of the largest to smallest singular value). The three Hessian matrices 
'Hrr, "^^cc ^cc essentially equivalent for investigating numerical issues and for testing 
whether a proposed minimizer of the second order expansion of /(r) = /(c) is a local (or even 
global) minimum. Thus, one can choose to work with the Hessian matrix which is easiest to 
compute and analyze. This is usually the c-complex Hessian H'^^, and it is often most convenient to 
determine numerical stability and optimality using H^^ even when the algorithm is being developed 
from one of the alternative perspectives (i.e., the real r or the c-real second order expansion). 

Now note that from (101) we immediately and easily have 

1 Ac^ Ac = 1 Ac^ S Ac = i {SAcf H^, Ac = i JKcf Ac = i Ac^ H^, Ac 

showing the equivalence of the c-real and c-complex second order terms in the expansion of /(c).^^ 
Combining this result with (98), we have shown the following equivalences between the second 
order terms in the various expansions of / under consideration in this note: 



2nd-Order Terms: ^ Ar^ Hrr ^^=\ ^cc = ^ Ac^ U^^ Ac = Re { Az^W^^ Az + Az^Ti;^^ Az} 



(104) 

where the second order expansion in r is given by (73), the c-complex expansion by (92), the 
expansion in terms of z by (99), and the c-real expansion by 



df(c) 1 

2nd-0rder Expansion in c-Real: /(c + Ac) = /(c) + ' Ac + - Ac"^ 'HlJc) Ac + h.o.t. 

oc 2 



(105) 

Note that all of the terms in (104) and (105) are real valued. 



''^One can show that the term Ac"^ TC^^ Ac is admissible if and only if 7^*^. = SM for M e £(C, C), which is the 
case here. 



K. Kreutz-Delgado — Copyright © 2003-2009, All Rights Reserved - Version UCSD-ECE275CG-S2009vl.O 40 



The expansion in of /(c) in terms of c-complex shown in (92) is not differentiable with respect 
to c as differentiation with respect to c-complex is not defined. (Recall, though, that we can differ- 
entiate the c-real expansion with respect to c-real.) However, (92) is differentiable with respect to 
z and z and can be viewed as a short-hand equivalent to the full (z, z) expansion provided by (99). 
Therefore, it is Equation (99) which is the natural form for optimization with respect to c-complex 
via a derivative-based approach, because only differentiation with respect to the components (z, z) 
of c-complex is well-posed. On the other hand, differentiation with respect to c-real is well-posed, 
so that one can optimize (105) by taking derivatives of (105) with respect to c-real itself. 

Note that (73), (92), and (105) are the natural forms to use for optimization via "completing 
the square." This is because the expansions in terms of r, c-complex, and c-real are less awkward 
for completing-the-square purposes than the expansion in z provided by (99).^^ Note, further that 
the expansions (73) and (92) both have a form amenable to optimization by completing the square 
and both are differentiable with respect to the expansion variable itself. 

The various second order expansions developed above can be found in references [36], [25] 
and [32]. In [25], van den Bos shows the equality of the first, second, and third second-order terms 
shown in equation (98) but does not mention the fourth (which, anyway, naturally follows from 
the third term in (98) via a simple further expansion in terms of z and z). Indeed, the approach 
used in this note is a more detailed elaboration of the derivations presented by van den Bos in 
[25]. In reference [32] Yan and Fan show the equality of the first and last terms in (98), but, while 
they cite the results of van den Bos [25] regarding the middle terms in (98), do not appear to have 
appreciated that the fourth term in (98) is an immediate consequence of the second or third terms, 
and instead derived it from scratch using an alternative, "brute force" approach. 

Quadratic Minimization and the Newton Algorithm. The Newton algorithm for minimizing a 
scalar function f{z) exploits the fact that it is generally straightforward to minimize the quadratic 
approximations provided by second order expansions such as (73), (92), (99), and (105). The 
Newton method starts with an initial estimate of the optimal solution, say c, then expands /(c) 
about the estimate c to second order in Ac = c — c, and then minimizes the resulting second 
order approximation of /(c) with respect to Ac. Having determined an estimated update Ac in 
this manner, one updates the original estimate c ^ c + aAc, for some small "stepsize" a > 0, 
and then starts the optimization cycle all over again. For appropriate choices of the stepsize a, this 
iterative approximate quadratic optimization algorithm can result in a sequence of estimates Cq, Ci, 
C2, ■ ■ • , which converges to the true optimal solution extremely quickly [34]. 

Note that the optimal solution to the quadratic approximations provided by (73), (92), and (105) 
can be immediately written down using the "completing-the-square" procedure assuming that the 
relevant Hessians are all invertible: 




(from the r expansion (73)) 



(106) 



Although (99) can also be optimized by completing the square. 
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in 



cc. 



-(Kc) 



V 9c 
dc 



H 



(from the c-complex expansion (92)) (107) 



(from the c-real expansion (105)) . 



(108) 



Solutions (106) and (107) can also be found in van den Bos [25]. Note that Ac*^ is an admissible 
solution, i.e., that 

Acc e C 

as required for self-consistency of our theory, as a consequence of the fact that and 
(H^J-^ satisfy 



(^)%C and (7^y-^G£(C,C) 



with the latter condition a consequence of property (82) and the fact that Ti^^. G C{C, C). If this 
were not the case, then we generally would have the meaningless answer that Ac*^ ^ C. 

The admissibility of the solution (108) follows from the admissibility of (107). This will be 
evident from the fact, as we shall show, that all of the solutions (106)-(108) must all correspond to 
the same update, 

Ac^= Ac« = JAr . 



Note that 



Ac^ 



Ifdfjc) 
dc 



H 



— I — JT^rrJ 



H 



2 dr 



H 



- (JKJ-T'J (^) 



T 



(from (67) and (90)) 
(from (63)) 



= -J(Kr) 
= JAr 
as required. On the other hand, 

Ac* = - 



I dr 



u \-lfdf{c) 
\ dc 



-(n 
in 



dc 



c \-l fdf jc) 
ccJ 



(from (101)) 



\ 9c 
c \-lfdfic) 



S 



V 9c 
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Thus, the updates (106)-(108) are indeed equivalent. 

The updates (106) and (108), determined via a completing the square argument, can alterna- 
tively be obtained by setting the (real) derivatives of their respective quadratically-approximated 
loss functions to zero, and solving the necessary condition for an optimum. Note that if we attempt 
to (erroneously) take the (complex) derivative of (92) with respect to c-complex and then set this 
expression to zero, the resulting "solution" will be off by a factor of two.^^ In the latter case, we 
must instead take the derivatives of (99) with respect to z and z and set the resulting expressions 
to zero in order to obtain the optimal solution. 

At convergence, the Newton algorithm will produce a solution to the necessary first-order con- 
dition 

^ = 0, 
dc 

and this point will be a local minimum of /(■) if the Hessians are strictly positive definite at this 
point. Typically, one would verify positive definiteness of the c-complex Hessian at the solution 
point c. 

As done in [36] and [32], the solution to the quadratic minimization problem provided by (106)- 
(108) can be expressed in a closed form expression which directly produces the solution z G C". 
To do so, we rewrite the solution (107) for the Newton update Ac as 



ULAc 



df{c) 



H 



which we then write in expanded form in terms of z and z 



zz 
zz 





(109) 



Assuming that Ti^^ is positive definite, then Ti^^ is invertible and the second block row in (109) 
results in 



Az = -H-jH.^/\7. - H-J ( ^ 



dz 



*^In a numerical solution procedure a constant factor error in the updates can be absorbed into the update step- 
size factor and therefore will likely not be noticed in simulations or applications. However, the claim that a specific 
step-size values results in stable or unstable convergence might not be confirmed in an experiment using the correctly 
computed updates. 

^^This is the procedure used in [36] and [32]. 
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Plugging this into the first block row of (109) then yields the Newton algorithm update equation 
where 

/Tzz '^-ZZ '^-ZZ '''ZZ '^zz 

is the Schur complement of TYzz in Ti^c- Equation (110) is equivalent to the solution given as 
Equation (A. 12) in [36]. Invertibility of the Schur complement Tizz follows from our assumption 
that is positive definite, and the Newton update is therefore given by 

^^ = [nzz-H-zzn-Jnz^~'i^H-zzH-j{^^^ (iii) 

The matrices Ti^z and Ti^x = {H7.7. — HzzH^iHzz,) in (110) are invertible if and only if TY^c is 
invertible. Note that invertibility of Hzz (equivalently, Hzz = 'Hzz) is not a sufficient condition for 
the Schur complement to be nonsingular. However, if Tizz = ^zz = then invertibility of Hzz is 
a necessary and sufficient condition for a solution Az to exist. 

As noted by Yan & Fan [32], the need to guarantee positive definiteness of the Schur comple- 
ment Hzz = {Hzz — HzzH^iHzz) is a significant computational burden for an on-line adaptive 
filtering algorithm to bear. For this reason, to improve the numerical robustness of the Newton 
algorithm and to provide a substantial simplification, they suggest making the approximation that 
the block off-diagonal elements of H'^^ are zero 

H^zz "H-zz ~ 
which results in the simpler approximate solution 

The argument given by Yan and Fan supporting the use of the approximation Hzz ~ is that as the 
Newton algorithm converges to the optimal solution z = zq, setting Hzz "to zero implies that we 
will use a quadratic function to approximate the cost near zq" [32]. However Yan and Fan do not 
give a formal definition of a "quadratic function" and this statement is not generally true as there 
is no a priori reason why the off-diagonal block matrix elements of the Newton Hessian should be 
zero, or approach zero, as we demonstrate in Example 2 of the Applications section below. 

However, as we shall discuss later below, setting the block off-diagonal elements to zero is 
justifiable, but not necessarily as an approximation to the Newton algorithm. Setting the block 
off-diagonal elements in the Newton Hessian to zero, results in an alternative, "quasi-Newton" 
algorithm which can be studied in its own right as a competitor algorithm to the Newton algorithm, 
the Gauss-Newton algorithm, or the gradient descent algorithm. 

^^That is not to say that there can't be conditions under which the quasi-Newton algorithm does converge to the 
Newton algorithm. Just as one can give conditions for which the Gauss-Newton algorithm converges to the Newton 
algorithm, one should be able to do the same for the quasi-Newton algorithm. 
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Nonlinear Least-Squares: Gauss vs. Newton. In this section we are interested in finding an 
approximate solution, z, to the nonlinear inverse problem 

g(z) ^ y 

for known y E C"^ and known real-analytic function g : C" — C™. We desire a least-squares 
solution, which is a solution that minimizes the weighted least-squares loss function^" 

^(z) = l\\y- g(z) \\w = liy- S{z)f W{y- g(z)) 

where 14^ is a Hermitian positive-definite weighting matrix. Although the nonlinear function g is 
assumed to be real-analytic, in general it is assumed to be not holomorphic (i.e., g is not complex- 
analytic in z). 

In the subsequent development we will analyze the problem using the c-real perspective devel- 
oped in the preceding discussions. Thus, the loss function is assumed to be re-expressible in terms 
of c, 

£(c) = i||y-g(c)f^ = i (y-g(c))^W^(y-g(c)) . (113) 

Intermediate quantities produced from this perspective^' may have a different functional form than 
those produced purely within the z E Z perspective, but the end results will be the same. 

We will consider two iterative algorithms for minimizing the loss function (113): The Newton 

algorithm, discussed above, and the Gauss-Newton algorithm which is usually a somewhat simpler, 

yet related, method for iteratively finding a solution which minimizes a least-squares function of 
the form (113).^2 

As discussed earlier, the Newton method is based on an iterative quadratic expansion and min- 
imization of the loss function £(z) about a current solution estimation, z. Specifically the Newton 
method minimizes an approximation to i{c) = i{z) based on the second order expansion of ^(c) 
in Ac about a current solution estimate c = col(z, z), 

i{c + Ac) ^ £(Ac)^™'™ 

where we define the Newton approximate loss function, 

£(Ac)^="'°" = £(c) + Ac + ^Ac^ nl,{c) Ac. (114) 

^"The factor of i has been included for notational convenience in the ensuing derivations. If it is removed, some 
of the intermediate quantities derived subsequently (such as Hessians, etc.) will differ by a factor of 2, although the 
ultimate answer is independent of any overall constant factor of the loss function. If in your own problem solving 
ventures, your intermediate quantities appear to be off by a factor of 2 relative to the results given in this note, you 
should check whether your loss function does or does not have this factor 

^'Such as the Gauss-Newton Hessian to be discussed below. 

^^The Newton algorithm is a general method that can be used to minimize a variety of different loss functions, while 
the Gauss-Newton algorithm is a least-squares estimation method which is specific to the problem of minimizing the 
least-squares loss function (113). 
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- Newton 



Minimizing the Newton loss function £(Ac)'^^"'°" then results in a correction Ac which is then 

— — Newton 

used to update the estimate c ^ c + aAc for some stepsize a > 0. The algorithm then starts 
all over again. As mentioned above, a "completing-the- square" argument can be invoked to readily 
show that the correction which minimizes the quadratic Newton loss function is given by 



- Newton 



di{c) 



H 



Ac =-Kc(6r^^J (115) 

provided that the c-complex Hessian ?^cc(^) is invertible. Because it defines the second-order 
term in the Newton loss function and directly enters into the Newton correction, we will often 
refer to 7Y^p(c) as the Newton Hessian . If we block partition the Newton Hessian and solve for 

Newton 

the correction Az , we obtain the solution (111) which we earlier derived for a more general 
(possibly non-quadratic) loss function. 

We now determine the form of the cogradient of the least-squares loss function (1 13). This 
is done by utilizing the c-real perspective which allows us to take (real) cogradients with respect 
to c-real. First, however, it is convenient to define the compound Jacobian G{c) o/g(c) as 

G'(e)^^^(¥ ^) = (j^(c) j|(c))ec-x2r (116) 



Setting e = y — g(c), we have 



73 



dc 2 dc 



- e"W— e + - e^W^— e 
2 dc 2 dc 

2 dc 2 dc 



le^WG --e^W^ (^S 



2 2 \dc 

= --e'^WG --e^W^GS 
2 2 

or 

^ = -Le''WG --e^WGS. (117) 
dc 2 2 

This expression for |^ is admissible, as required, as it is readily verified that 



If)" -'if 

dc J \dc 



as per the requirement given in (76). 



73 



Remember that ^ is only well-defined as a derivative within the c-real framework. 
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The linear term in the Newton loss function i'^™'™ is therefore given by 

f)f 1 1 

—Ac = — e^WG Ac - -e^WGS Ac 
dc 2 2 

= --e^W^GAc- - e^WG Ac 
2 2 

= -Re {e^iyCAc}. 

Thus 

— Ac = -Re {e^VTGAc} = -Re {(y - g(c))^ W^G Acj . (118) 

If the reader has any doubts as to the validity or correctness of this derivation, she/he is invited to 
show that the right-hand side of (1 18) is equal to 2 Re | ^ Az} as required from equation (86). 

Before continuing on to determine the functional form of the c-complex Hessian 7-^^^, (c) needed 
to form the Newton loss function and solution, we turn first to a discussion of the Gauss-Newton 
algorithm. 

Whereas the Newton method is based on an iterative quadratic expansion and minimization of 
the loss function £(z) about a current solution estimation, z. The Gauss-Newton method is based 
on iterative "relinearization" of the system equations y ^ g(z) about the current estimate, z and 
minimization of the resulting approximate least-squares problem. We put "linearization " in quotes 
because (unless the function g happens to be holomorphic) generally we are not linearizing g with 
respect to z but, rather, we are linearizing with respect to c = co/(z, z). 

Expanding the system equations y ^ g(z) about a current estimate z, we have 
y - g(z) = y - g(z + Az) ^ y - (^g(z) + ^Az + ^Az^ 

where Az = z — z and Az = Az = z — z = z — z. Note that the approximation to g is not a linear 
function of z as complex conjugation is a nonlinear operation on z. However, if g is holomorphic, 
then ll = 0, in which case the approximation is linear in z. Although the approximation of g 
generally is not linear in z, it is linear in c = col(z, z), and we rewrite the approximation as 

y - g(c) = y - g(c + Ac) ^ Ay - G(c) Ac (1 19) 

where Ay = y — g(z), c = col(z, z), Ac = c — c, and G{c) is the (compound) Jacobian mapping 
of g evaluated at the current estimate c given in Equation (116). With this approximation, the loss 
function (1 13) is approximated by the following quadratic loss function (notationally suppressing 
the dependence on c), 

i{c) = i{c + Ac) ^ i{Acy'"'' 
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where 

i{Acr- = i||Ay-GAc||J^ 

= ^{Ay-GAcfWiAy-GAc) 

= -||Ayf -Re {Ay^W G Ac} + -Ac^ G^WG Ac 

= i(c) + Ac + -Ac^ G^WG Ac. (from (118) 
oc 2 

Unfortunately, the resulting quadratic form 

OC 2 

is no? admissible as it stands. This is because the matrix G^WG is not admissible, 
This can be seen by showing that the condition (78) is violated: 



SG^WGS 




Fortunately, we can rewrite the quadratic form (120) as an equivalent form which is admissible 
on C. To do this note that G^WG is Hermitian, so that 

Ac" G^WG Ac = Ac^G^WGAc e M . 

Also recall from Equation (79) that P{G"WG) G C{C,C) and Ac e C ^ SAc = Ac. For an 

^'^And thus the complex Gauss-Newton algorithm is generally more complicated in form than the real Gauss-Newton 
algorithm for which the quadratic form (120) is meaningful. 
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admissible variation Ac G C we have^^ 

Ac^G^WGAc = Ac^P{G^WG)Ac + Ac^ [G^WG - P{G^WG)) Ac 
= Ac^P{G^WG)Ac + ^Ac^ (g^WG - SG^WGS^ Ac 

= Ac^P{G^WG)Ac + i (^Ac^G^WGAc - Ac^G^WGAc 

= Ac^P{G^WG)Ac + 
= Ac^P{G^WG)Ac. 

Thus we have shown that on the space of admissible variations, Ac G C, the inadmissible 
quadratic form (120) is equivalent to the admissible quadratic form (the Gauss-Newton approxi- 
mate loss function) 

i{Acy- = £(c) + Ac + ^Ac^ Kric) Ac (121) 



where 



Gauss-Newton Hessian n°ric) = P {G^ic)WG{c)) 



(122) 



denotes the Gauss-Newton Hessian. Note that the Gauss-Newton Hessian is the exact Hessian 
matrix of the Gauss-Newton approximate loss function. 

Note that the Gauss-Newton Hessian 7Y°™''(c) is Hermitian and always guaranteed to be at least 
positive semi-definite, and guaranteed to be positive definite if g is assumed to be one-to-one (and 
thereby ensuring that the compound Jacobian matrix G has full column rank). This is in contrast 
to the Newton (i.e., the c-complex) Hessian 'H'^^{c) which, unfortunately, can be indefinite or rank 
deficient even though it is Hermitian and even if g is one-to-one. 

Assuming that '^^^''(c) is invertible, the correction which minimizes the Gauss-Newton ap- 
proximate loss function (121) is given by 

Ac""' = ~Kricr'(^Y . (123) 



dc 



Because of the admissibility of T^cc ' ^^'^ ' resulting solution is admissible Ac G C. 

Comparing Equations (115) and (123), it is evident that the difference between the two al- 
gorithms resides in the difference between the Newton Hessian, 7-^^j.(c), which is the actual c- 
complex Hessian of the least-squares loss function i{c), and the Gauss-Newton Hessian 7-^cc"'(^) 
which has an as yet unclear relationship to £(c).^^ For this reason, we now turn to a discussion of 
the relationship between the Hessians 'H'^^{c) and 7Y°^"''(c). 



H 



~ Gaus; 



^^Note that the ensuing derivation does not imply that G^WG = P{G^WG), a fact which would contradict our 
claim that G^WG is not admissible. This is because in the derivation we are not allowing arbitrary vectors in C^" 
but are only admitting vectors Ac constrained to lie in C, Ac g C C C^". 

^^Note that, by construction, 7Y°™''*(c) is the Hessian matrix of the Gauss-Newton approximate loss function. The 
question is: what is its relationship to the least-squares loss function or the Newton approximate loss function? 
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We can compute the Newton Hessian H^^. from the relationship (see Equation (101)) 

n_/C _ C-iyK _ c ^ ( 

where ^ is taken to be a c-real cogradient operator. Note from (1 17) that, 

V^c) =-lf''We--SGHWe = -{B + SB), (124) 

where 

B = -G"We (125) 

with e = y — g(c). This results in 



dc I \dc / 2 



{B + SB) 



Also note that 



We have 



or 



dB__ dB_ _ fdB_ g\ _ dB_ ^ 
dc dc \dc J dc 



d fdtV 1 f ^dB dB 



^"'^ dc[dc) 2V dc~^ dc 



^- = d-c[d-c =21^+^ ^ ^ 



This yields 



«L = s«L = i(f + sf SI (127) 



with B given by (125), which we can write as 

Kc = 5Kc = p(^) • (128) 

Recall that H"^^ must be admissible. The function P(-) produces admissible matrices which map 
from C to C, and thereby ensures that the right-hand side of equation (128) is indeed an admissible 
matrix, as required for self-consistency of our development. The presence of the operator P does 
not show up in the real case (which is the standard development given in textbooks) as ^ is 
automatically symmetric as required for admissibility in the real case. 

Note that B can be written as 
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where gi and [W (y — g) ],; denote the i-th (scalar) components of the vectors g and We = W{y — 

dB 
dc 



g) respectively. We can then compute ^ as 



oB / c'g \ I / a(7j 



or 



dB 
dc 



Equations (128) and (129) result in the following succinct relationship between the complex 
Newton and Gauss-Newton Hessians, 



Newton Hessian = = 7^°^"^^ - ^ 



i=l 



(130) 



where the Gauss-Newton Hessian 7^°™'' 1^ given by (122) and 



Equation (130), which is our final result for the structural form of the Newton Hessian H^^, looks 
very much like the equivalent result for the real case.''^ The first term on the right-hand-side of 
(130) is the Gauss-Newton Hessian 7i°^"'', which is admissible, Hermitian and at least positive 
semidefinite (under the standard assumption that W is Hermitian positive definite). Below, we 
will show that the matrices W^c^ i = 1, ■ ■ ■ ,m, are all individually admissible and Hermitian.^** 
While the Gauss-Newton Hessian is always positive semidefinite (and always positive definite if 
g is one-to-one), the presence of the second term on the right-hand-side of (130) can cause the 
Newton Hessian to become indefinite, or even negative definite. 

We can now understand the relationship between the Gauss-Newton method and the Newton 
method when applied to the problem of minimizing the least-squares loss function. The Gauss- 
Newton method is an approximation to the Newton method which arises from ignoring the second 
term on the right-hand-side of (130). This approximation is not only easier to implement, it will 
generally have superior numerical properties as a consequence of the definiteness of the Gauss- 
Newton Hessian. Indeed, if the mapping g is onto, via the Gauss-Newton algorithm one can 
produce a sequence of estimates c^, /c = 1, 2, 3, ■ ■ ■ , which drives e(ck) = y — g(cfc), and hence 



^'The primary difference is due to the presence of the projector P in the complex Newton algorithm. Despite the 
similarity, note that it takes much more work to rigorously derive the complex Newton- Algorithm! 

^^Of course, because TY^c ™d ^cc"'" ^re both Hermitian and admissible the total sum ^cc must be Hermitian 
and admissible. 
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(with some additional smoothness assumptions on g) the second term on the right-hand- side of 
(130), to zero as A; ^ oo. In which case, asymptotically there will be little difference in the 
convergence properties between the Newton and Gauss-Newton methods. This property is well 
known in the classical optimization literature, which suggests that by working within the c-real 
perspective, we may be able to utilize a variety of insights that have been developed for the Newton 
and Gauss-Newton methods when optimizing over real vector spaces. 

We will now demonstrate that each individual term 7i^^, i = 1, ■ ■ ■ , m, in (130) is admissible 
and Hermitian. Note that the "raw" matrix 



[We_ 



d_ 
dc 



dc 



H 



is neither Hermitian nor admissible because of the presence of the complex scalar factor [W^e]^. 
Fortunately, the processing of the second matrix of partial derivatives by the operator P to form 



the matrix W^l via 



' '-cc 



P [Wi 



d_ 
dc 



dc 



creates a matrix which is both admissible and Hermitian. The fact that is admissible is obvious, 
as the projector P is idempotent. We will now prove that H 

Define the matrix 



is Hermitian. 



^cc (fi'i 



d_ 
dc 



dc 



H 



(132) 



and note that 



d_ 
dc 



dgi 
dc 



H' 



H 



d_ 
dc 



dgi 
dc 



d_ 
dE 



dgi 
dc 



d_ 
dc 



dgi 
dc 



H 



which shows that Acc{gi) has the property that 



(133) 



Now note that 



d f dg 



dc \ dc 



H 



d f dgi 



dc V dc 



H 



d_ 
dc 



S 



d ( dg, 



dgi 
dc 



H' 



H 



dc \ dc 

d_ f dgi' " 
dc \ dc 
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which establishes the second property that 

SA^MS = Aec(^?.) . (134) 
Finally note that properties (133) and (134) together yield the property 



Setting tti = [We]., we have 



n,(i)_-n^ A I ^^ _ ai Accigi) + S Accigi) S _ ai Accigi) + ai S Accigi) S _ ai A^ciQi) + ai Accigi)" 

'Icc — -t^l^i ^cc[gi)) — ^ — ^ — ^ 

which is obviously Hermitian. Note that the action of the projector P on "raw" matrix Acdgi), 
Hermitian symmetrizes the matrix Acdgi)- 

Below, we will examine the least- squares algorithms at the block-component level, and will 
show that significant simplifications occur when g(z) is holomorphic. 

Generalized Gradient Descent Algorithms. As in the real case, the Newton and Gauss-Newton 
algorithms can be viewed as special instances of a family of generalized gradient descent algo- 
rithms. Given a general real- valued loss function £(c) which we wish to minimize^^ and a current 
estimate, c of optimal solution, we can determine an update of our estimate to a new value c„,„ 
which will decrease the loss function as follows. 

For the loss function £(c), with c = c + dc, we have 

di(c) = i(c + dc) - i(c) = ^^dc 

oc 

which is just the differential limit of the first order expansion 

Ai(c; a) = i(c + aAc) - i(c) ^ a^^Ac . 

oc 

The stepsize « > is a control parameter which regulates the accuracy of the first order approxi- 
mation assuming that 

a — ^ ^ aAc — * dc and A£(c; a) — > di{c) . 

If we assume that C is a Cartesian space,^° then the gradient of i{c) is given by^^ 

'd£{c)^ ^ 



Ve^(c) 



dc 



^'^The loss function does not have to be restricted to the least-squares loss considered above. 

'^'^I.e., We assume that C has identity metric tensor. We call the resulting gradient a Cartesian gradient (if the metric 
tensor assumption fl^ = I is true for the space of interest) or a naive gradient (if the identity metric tensor assumption 
is false, but made anyway for convenience). 

**'Note for future reference that the gradient has been specifically computed in Equation (124) for the special case 
when ^(c) is the least-squares loss function (113). 
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Take the update to be the generalized gradient descent correction 

Ac = -Q{c) = -g(c) Vc£(c) (135) 

where Q{c) is a Hermitian matrix function of c which is assumed to be positive definite when 
evaluated at the value c.^^ This then yields the key stability condition^^ 

A£(c;a) ^ -a\\VJ{c)\\l ^ -aVJ{cfQVJic) < 0, (136) 

where the right-hand-side is equal to zero if and only if 

Vc^(c) = . 



Thus if the stepsize parameter a is chosen small enough, making the update 

Cnew = c + a Ac = c-Q Vc^(c) 

results in 

£(c_) = £(c + aAc) = £(c) + A£(c; a) ^ i{c) - a|| Ve^(c) ||^ < £(c) 

showing that we either have a nontrivial update of the value of c which results in a strict decrease 
in the value of the loss function, or we have no update of c nor decrease of the loss function 
because c is a stationary point. If the loss function £{c) is bounded from below, iterating on this 
procedure starting from a estimate Ci will produce a sequence of estimates Ci, i = 1, 2, 3, ■ ■ ■ , 
which will converge to a local minimum of the loss function. This simple procedure is the basis 
for all generalized gradient descent algorithms. 

Assuming that we begin with an admissible estimate, Ci, for this procedure to be valid, we 
require that the sequence of estimates Cj, z = 1, 2, 3, ■ ■ ■ , be admissible, which is true if the corre- 
sponding updates Ac are admissible, 

Ac = -g(cO Vc,^(cO = -Q(cO eC, z = 1, 2, ■ ■ ■ . 

We have established the admissibility of Vc^(c) = G C above. It is evident that in order 

for a generalized gradient descent algorithm ( GDA ) to be admissible it must be the case that Q be 
admissible, 

Generalized GDA is Admissible <^ Generalized Gradient Q-Matrix is Admissible, Q G C{C,C) . 

^^The fact that Q is otherwise arbitrary (except for the admissibihty criterion discussed below) is what makes the 
resulting algorithm a generalized gradient descent algorithm. When Q = I, we obtain the standard (naive) gradient 
descent algorithm. 

'^^We interpret the stability condition to mean that for a small enough stepsize a > 0, we will have A£{c; a) < 0. 
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Furthermore, a sufficient condition that the resulting algorithm be stablizable^^ is that Q be Her- 
mitian and positive definite. Note that given a candidate Hermitian positive definite matrix, Q' , 
which is not admissible, 

we can transform it into an admissible Hermitian positive definite matrix via the projection 

Q = V{Q')eC{C,C). 
It can be much trickier to ensure that Q remains positive definite under the action of P. 
If we set 

/nNewton/'^\ ri/ Newton /' ^ M — 1 

with 

then we obtain the Newton algorithm (115). If we take the loss function to be the least-squares loss 
function (113) and set 

we obtain the Gauss-Newton algorithm (123). Whereas the Gauss-Newton algorithm generally 
has a positive definite Q-matrix (assuming that g{c) is one-to-one), the Newton algorithm can 
have convergence problems due to the Newton Hessian TY^™'"" = T^cc becoming indefinite. Note 
that taking 

Q = l. 

which we refer to as the "Cartesian," "standard," "simple," or "naive" choice (depending on the 
context) results in the standard gradient descent algorithm which is stable for a small enough 
stepsize so that the stability condition (136) holds. 

An important practical issue is the problem of stability versus speed of convergence. It is well- 
known that the Newton algorithm tends to have a very fast rate of convergence, but at the cost of 
constructing and inverting the Newton Hessian 7-^^™'°" = TY^c ^iid potentially encountering more 
difficult algorithm instability problems. On the other hand, standard gradient descent {Q = I) 
tends to be very stable and much cheaper to implement, but can have very long convergence times. 

The Gauss-Newton algorithm, which is an option available when the loss function £(c) is the 
least-squares loss function (113), is considered an excellent trade-off between the Newton algo- 
rithm and standard gradient descent. The Gauss-Newton Hessian 7^°^ generally simpler in form 
and, if g(c) is one-to-one, is always positive definite. Furthermore, if g(c) is also onto, assuming 
the algorithm converges, the Gauss-Newton and Newton algorithms are asymptotically equivalent. 

We can also begin to gain some insight into the proposal by Yan and Fan [32] to ignore the 
block off-diagonal elements of the Newton Hessian,^^ 

z 



/Newton '1-/^ 

'-CC '''CC 



z 



I.e., that a small enough step size can be chosen to ensure that the stability condition (136) is satisfied. 
**^The values of the block elements of 7Y^c"'°" will be computed for the special case of the least-squares loss function 
(113) later below. 
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As mentioned earlier, Yan and Fan make the claim in [32] that the block off-diagonal elements 
vanish for a quadratic loss function. As noted above, and shown in an example below, this is 
not generally true.^^ However, it is reasonable to ask what harm (if any), or what benefit (if any) 
can accrue by constructing a new^^ generalized gradient descent algorithm as a modification to 
the Newton algorithm created by simply ignoring the block off-diagonal elements in the Newton 
Hessian and working instead with the simplified quasi-Newton Hessian, 



quasi-Newton _A_ njC ^ { ^ZZ 



1_/quasi-Newton ^ n / 



n 



zz 



This results in a new generalized gradient descent algorithm, which we call the quasi-Newton 
algorithm, which is somewhere in complexity between the Newton algorithm and standard gradient 
descent. Note that the hermitian matrix Tizz is positive definite if and only if Tizz is positive 
definite. Thus invertibility and positive-definiteness of the quasi-Newton Hessian 7-^^"^" = 7i^^ 
is equivalent to invertibility and positive definiteness of the block element Tizz- 

On the other hand, invertibility and positive definiteness of TYzz is only a necessary condition 
forinvertibility and positive definiteness of the complete Newton Hessian 7i^^"'°" = Ticc- Assuming 
that is positive definite, we have the well-known factorization 

I o\.,, // -n^,n-:\ /^zz 



-K.Hz-i i) ^'^'^ / 'J - V nj ^^^^^ 

where 

'T-zz 'T-zz ''-zz''-zz ''-zz 

is the Schur complement of TYzz in ^cc- From the factorization (137) we immediately obtain the 
useful condition 

rank (H^J = rank {H^^) + rank {u^^ . (138) 

Note from condition (138) that the Newton Hessian 7-^^™'™ = 7-^^^ is positive definite if and 
only if TYzz <^nd its Schur complement TYzz are both positive definite. Thus it is obviously a more 
difficult matter to ascertain and ensure the stability of the Newton Hessian than to do the same for 
the quasi-Newton Hessian. 



The quasi-Newton algorithm is constructed by forming the Q matrix from the quasi-Newton 

^zz' 



Hessian 7^^."-' n™'™ = 



Q Pseudo-Newton /'T/qiiasi-NewtoiiN ~1 / nj'^ 



which is admissible and hermitian, and positive definite provided TYzz = ^zz is positive definite. 
Thus, if TYzz = ^zz is positive definite, the quasi-Newton algorithm is guaranteed to be stable 



What is true, as we've noted, is that for a quadratic loss function, the Gauss-Newton and Newton Hessians asymp- 
totically become equal. 

'^^I.e., no approximation algorithms are invoked. 
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(assuming a small enough stepsize a > so that the stability condition (136) is satisfied). With 
this choice of Q in (135), the quasi-Newton update is given by^*^ 

which is just the simplification shown earlier in Equation (112) and proposed by Yan and Fan in 
[32]. However, unlike Yan and Fan, we do not present the quasi-Newton algorithm as an approx- 
imation to the Newton algorithm, but rather as one more algorithm in the family of generalized 
Newton algorithms indexed by the choice of the matrix Q. 

Indeed, recognizing that the Gauss-Newton algorithm potentially has better stability properties 
than the Newton algorithm, naturally leads us to propose a quasi-Gauss -Newton algorithm for 
minimizing the least-squares lose function (113) as follows. Because the hermitian Gauss-Newton 
Hessian is admissible, it can be partitioned as 

-^Gauss _ ( ^ZZ ^ZZ 



with f/gz = f/^.^^ The Gauss-Newton Hessian is positive-definite if and only if U^^ (equivalently 
t/zz) and its Schur complement U^t. = U^z — U^xU^z U^z are invertible. 

On the other hand the quasi-Gauss-Newton Hessian, 

^quasi-Gauss A / ^ZZ 



' Uzz, 

is positive definite if and only if Uzz is positive definite. Choosing 

jj-^ 

quasi-Gauss / 'T/quasi -Gauss N ~ 1 j 



''CC 



^^quasi-Gauss 

results in the quasi-Gauss-Newton algorithm 



f/.. ' 



^^jq-i-G-- = _f/-i [ ^ ) (140) 



7Z 



H 



which is guaranteed to be stable (for a small enough stepsize so that the stability condition (136) 
is satisfied) if Uzz is positive definite. 

Note that TYzz can become indefinite even while Uzz remains positive definite. Thus, the quasi- 
Gauss-Newton algorithm appears to be generally easier to stabilize than the quasi-Newton algo- 
rithm. Furthermore, if g is onto, we expect that asymptotically the quasi-Gauss-Newton and quasi- 
Newton algorithm become equivalent. Thus the quasi-Gauss-Newton algorithm is seen to stand in 



**^We can ignore the remaining update equation as it is just the complex conjugate of the shown update equation. 
'^''The values of these block components will be computed below. 
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the same relationship to the quasi-Newton algorithm as the Gauss-Newton algorithm does to the 
Newton algorithm. 

Without too much effort, we can construct the block matrix components needed to implement 
the Newton and Gauss-Newton algorithms developed above in order to minimize the least-squares 
loss function (113).^° 

Let us first look at the elements needed to implement the Gauss-Newton algorithm. From 
Equation (122) and the derivations following Equation (120) one obtains 




which is positive definite, assuming that W is positive definite and that g is one-to-one. Similarly, 

one finds that 

dg\ f dg 




Also Uzz = Uzz and U^z = Uzz- We have now completely specified the Gauss-Newton Hessian 
-^Gauss quasi-Gauss-Newton Hessian at the block components level, 

T/Gauss I Uzz Uzz \ T/quasi-Gauss A ( U ZZ 

~ \ Uz^ U-zz) " V Uzz 

Now note the important fact that Uzz = Uzz = when g is holomorphicl Thus, when g is 
holomorphic there is no difference between the Gauss-Newton and pseudo-Gauss-Newton algo- 
rithms.^' Furthermore, when g(z) is holomorphic, Uzz simplifies to 

where Jg is the Jacobian matrix of g. 

Now let us turn to the issue of computing the elements need to implement the Newton Algo- 
rithm, recalling that the Newton Hessian is block partitioned as 



'T/Newton n_l' 

''-CC '''CC 



C / ^zz T~^zz 

'Hzz 'H zz 



One can readily relate the block components Hzz and Hzz to the matrices Uzz and Uzz used in the 
Gauss-Newton and quasi-Gauss-Newton algorithms by use of Equation (130). We find that 



"^^zz = Uzz — ^ K 



ZZ 



i=l 



'"This, of course, results in only a special case application of the Newton and quasi-Newton algorithms, both of 
which can be applied to more general loss functions. 

^'Recall that g(z) is holomorphic (analytic in z) if and only if the Cauchy-Riemann condition = is satisfied. 
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with 



y(') = - 
2 



d [ dgi{z) 



H 



d f %(z) 



9z \ 9z 

where e = y — g(z). Similarly, we find that 
and 



dz \ dz 



H 



(144) 



(i) 
ZZ 



i=l 



r(i) 
' ZZ 



1 

2 



d ( %(z) 



dz \ dz 



H 



d f %(z) 



dz \ dz 



H 



[We] 



(145) 



Furthermore, Vgz = Kz and V^z = Vg^. 



Note that neither V^z nor V^z vanish when g is holomorphic, but instead simplify to 



1 d_ 

2 di 



%(z) 
dz 



H 



[We] 



(0 



and Viz 



1 9_ 

2 ^ 



%(z) 
dz 



H 



[We] 



(146) 



We have shown that the relationship between the Newton Hessian and Gauss-Newton Hessian 
is given by 

T~^ZZ T~^ZZ \ ( ^ZZ Uz2 

. T~(-zz T~(-zz J \ ^zz Uzz , 



' zz 



0') 

' ZZ 



V'^' V- 



-^Newton -^Gauss 

In the special case when g(z) is holomorphic, the relationship becomes 

/ d (dg,{zr" 



f 'Hzz 


Ti-zz^ 


Uzz 


M 


\yizz 


'Hzz ) 


V 


f/zj 



■E 



dz V dz 



[We] 



oz \ dz I 



9 /%i(z) 



-^Newton 



\dz 



dz 



[We] 



d_ / dg,(z) 
dz V dz 



[We] 



This shows that if g(z) is holomorphic, so that the block off-diagonal elements of the Gauss- 
Newton Hessian vanish, and g(z) is also onto, so that asymptotically we expect that e 0, then 
the claim of Yan and Fan in [32] that setting the block off-diagonal elements of the Hessian matrix 
can proved a good approximation to the Hessian matrix is reasonable, at least when optimizing the 
least-squares loss function. However, when e ^ the Newton least-squares loss function (114) 
reduces to the Gauss-Newton loss function (121), so that in the least-squares case one may as 
well make the move immediately to the even simpler Gauss-Newton algorithm (which in this case 
coincides with the quasi-Gauss-Newton algorithm). 

However, the real point to be made is that any generalized gradient descent algorithm is worthy 
of consideration,^^ provided that it is admissible, provably stable, and (at least locally) convergent 



92 



I.e., we don't have to necessarily invoke an approximation argument. 
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to the desired optimal solution. After all the standard gradient descent algorithm corresponds to 
the cheapest "approximation" of all, namely that 

•-^.^Newton J 

and very few will deny the utility of this algorithm, even though as an "approximation" to the 
Newton algorithm it might be far from correct. The resulting algorithm has intrinsic merit as an 
algorithm in its own right, namely as the member of the family of gradient descent algorithms 
corresponding to the simplest choice of the Q-matrix, 

Q = I. 

In the end, if the algorithm works, it's ok. As it is said, "the proof is in the pudding."^^ 

We see, then, that we have a variety of algorithms at hand which fit within the framework of 
generalized gradient descent algorithms. These algorithms are characterized by the specific choice 
of the Q-matrix in the gradient descent algorithm, and include (roughly in the expected order 
of decreasing complexity, decreasing ideal performance, and increasing stability when applied to 
the least-squares loss function): 1) the Newton algorithm, 2) the quasi-Newton algorithm, 3) the 
Gauss-Newton algorithm, 4) the quasi-Gauss-Newton algorithm, and 5) standard gradient descent. 
Note that the Newton, quasi-Newton, and standard gradient descent algorithms are algorithms for 
minimizing a general loss function, while the Gauss-Newton and quasi-Gauss-Newton algorithms 
are methods for minimizing the least-squares loss function (113). 

For convenience, we will now summarize the generalized gradient descent algorithms that we 
have developed in this note. In all of the algorithms, the update step is given by 

c ^ c + aAc 



or, equivalently, 

z ^ z + aAz 

for a specific choice of the stepsize « > 0. The stability claims made are based on the assumption 
that a has been chosen small enough to ensure that the stability condition (136) is valid. Further- 
more, we use the shorthand notation 

and 

e(c) = y - g(c) . 

Note that in the taxonomy given below only the Newton Hessian Ti^J."'™ is generally the true Hes- 
sian of the loss function. 



Of course, we are allowed to ask what the performance of the Q = I standard gradient-descent algorithm is 
relative to the QNe»ton algorithm. 
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1 . Standard (a.k.a. Simple, Cartesian, or Naive) Gradient Descent. 

Applies to any smooth loss function which is bounded from below. 



^standar 



-V,^(c) 



dc 
dz 



H 



H 



Application to Least-Squares Loss Function (113): 

(f = -iG^'We - IS'mV^ = \ {B{c) + SBic)) 
where B{c) = -G{c)"We{c) 

y^ ^standard 



(f) 

y..standLU-d 

g(z) holomorphic: 



5(c) + SB{c] 

(z) 



^) Weiz) + {^] Weiz 



dz 



H 



H 



#) Weiz) 



H 



^Yweiz) 



We{z) 
ite)Ve(z) 



(f) 

y,standai-d 

Generally stable but slow. 

2. Gauss-Newton Algorithm. 

Applies to the least-squares loss function (113). 



Uzz Uzz 
Uzz Uzz 



where Uzz is given by (141), Uzz = Uzz, Uzz is given by (142), and Uzz = Uz 



dc 



H 



where 



(f )" = -iG^W^e - \SG^We = \ {B{c) + SB{c)) 
with B{c) = -G{c)^We{c) 

Az- = {Uzz - U-zzU,iUz^-' {u-zzU,^ (i)^ - (f )^ 



where 



(f) 



H 



dz 



H 



4^) Weiz) 



( dV\^ — (dn 
\dz) Kdz) 



H 
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g(z) holomorphic: 

takes the simpler form (143), f/^z = t4z, and Uzz = U^z 

u,,) 

H 



0. 



H 



=Uzz\^r= 



Az 



dz 



H 



g(21 



H 



We(z) 



Stability generally requires positive definiteness of both U^z and its Schur complement: 
Uzz = Uzz — UzzUz^Uzz- The need to step for positive-definiteness of the Schur complement 
can significantly increase the complexity of an on-line adaptive filtering algorithm. 

If g(z) is holomorphic, then stability only requires positive definiteness of the matrix Uzz = 
(^) ^ (%?^) ' which will be the case if g(z) is one-to-one. Thus, the algorithm may 
be easier to stabilize when g(z) is holomorphic. 
Convergence tends to be fast. 

3. Pseudo-Gauss-Newton Algorithm. 

Applies to the least-squares loss function (113). 



zz 
Uzz 



where Uzz is given by (141) and Uzz = Uz 



pseudo-Gauss / ^\ 



Q 



^pseudo-( 



[H 



pseudo-Gauss z' 

cc V^)\ 



U: 



Uzz 



-Q 



pseudo-Gausj 



de(c) 



H 



where 



(i)^ = -^G^W^e - ^SG^We = i (5(c) + SB{c)) 
with B{c) = -G{c)^We{c) 



Az' 



pseudo-Gauss 



[Uzzm 



-1 ( dejz) 

dz 



H 



n -1 



Where 



H 



We(z) 



(de\ 

\dz> 2 \ dz j " ' V dz 

g(z) holomorphic: 

Uzz takes the simpler form of (143) , and Uzz = Uzz- 

-iypsei,do-Gauss/A\ ^ ( ^ \ = ^ { ^ 
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(ir=-Kf^)"w^e(z 

y,pseudo-Gauss 



Stability requires positive definiteness of t/zz(z) = W j which will be the 

case if g(z) is one-to-one. 

Convergence is expected to be quick but generally slower than for Gauss-Newton due to loss 
of efficiency due to neglecting the block off-diagonal terms in the Gauss-Newton Hessian 
(off-set, however, by reduced complexity and possible gains in stability), except for the case 
when g(z) is holomorphic, in which case the two algorithms coincide. 

4. Newton-Algorithm. 

Applies to any smooth loss function which is bounded from below. 

'iyNewton|'A\ _ | '^Zz(c) "HzZ (c) 

Az-"- = (Kz - n,.n^iH.,y' {H,.H^i{§)" - mf} 

Application to the Least-Squares Loss Function (113): 

" V^zz -Hzzj'Kuzz ^'='\y^ 

_ -X7Gauss/;^\ V^m / I'zz V'zz 
— rice \^) - 2^i=l I yii) Y_{i) 

\ "^ZZ "^ZZ / 

f/zz is given by (141), U^z = Uzz, Uzz is given by (142), Uzz = 



is given by (144), V^^ = V^t is given by (145), K'l = 

H 



(f = -iG^'We - \SG^We = \ {B{c) + 55(c)) 
with B{c) = -G{c)^We{c) 

Az-"- = {Hzz - n-zzHiinz^y' {nzznii{§f - where 



(f )" = -I [(^?)"w^e(z) + (#) Ve(z) J; (§)" = (f ) 
g(z) holomorphic: 

'lyNewton _ I '-^zz >J \ v^»7i / ''zz "^zz \ _ n_/pxudo-aauss ( A\ sr^rn I ^ zz ^ zz 



H 
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zz 



(i) 



V^'J and take the simpler forms of (146), 
Uzz takes the simpler form of (143), f/gz = f4z 
Az^-™ = {Hzz - HzzH^inz^y' {nzzn^i{§f - {§f] where 



(i) 
ZZ 



(f) 



H 



m. 

dz 



H 



Stability generally requires positive definiteness of both Hzz and its Schur complement 
'^^zz = ('^^zz — 'Hzz'H^i'Hzz) ■ The need to step for positive-definiteness of the Schur com- 
plement can significantly increase the complexity of an on-line adaptive filtering algorithm. 

When minimizing the least-squares loss function, we expect stability to be greater when 
g(c) is holomorphic. This is particularly true if g(c) is also onto and the algorithm is con- 
vergent, as we then expect the difference between the Newton and Gauss-Newton Hessians 
(and hence the difference between the Newton and Gauss-Newton algorithms) to become 
negligible asymptotically. 

The Newton algorithm is known to have very fast convergence properties, provided it can be 
stabilized. 



5. Pseudo-Newton Algorithm. 

Applies to any smooth loss function which is bounded from below. 



•^.^pseudo-r 



zz V 









pseudo-Newton ^ ^ '' 
^psedudo-Newton 

^^2jpseudo-Newton _ 



pseudo-Newton /' A^l ^ 
CC \^)\ 



[K 

^^pseudo-Newton 



dijc) 
dc 



H 



[^zz(z)]-^ » 



H 



Application to the Least-Squares Loss Function (113): 



•^.^pseudo-Newton 



nzz{c) 

H,g(c) 



V 



1=1 





(i) 
ZZ 



-^pseudo-G.u,,s(^g^ 



T.Vz 

i=l 



(0 
zz 



V 

V^^ is given by (144) and 









m 
i=l 



Uz 





m 

Z-j ' zz 

i=l 



y^ ^pseudo-f 

with Bic 



-Q 



pseudo-Newton / 



Kz . Uzz is given by (141) and Uzz 

H 

where 



Uz 



dijc) 
dc 



^SG^We = I {B{c) + SB{c)) 
G{c)"We{c) 
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y-pseudo-f 



[Hzz(z)]" 



(2) 



H 



zz 



8£(z) 

az 



where 



f))''iye(z) + fe^'' 



g(z) holomorphic =^ 
f/zz takes the simpler form of (143), V-,, 
V^2 takes the simpler form (146), 



(i) 
ZZ 



Wz' 



a£(z) 

5z 



m 

zz I'zz 



H 



j=i 



lye(z) 

H 



dz 



[We] 



dz 



H 



We{z) 



Stability generally requires positive definiteness of Hzz- 

The pseudo-Newton is expected to be fast, but have a loss of efficiency relative to the Newton 
algorithm. When g(z) is holomorphic and onto, we expect good performance as asymptoti- 
cally a stabilized pseudo-Newton algorithm will coincide with the Newton algorithm. If g(z) 
is nonholomorphic, the pseudo-Newton and Newton algorithms will not coincide asymptot- 
ically, so the speed of the pseudo-Newton algorithm is expected to always lag the Newton 
algorithm. 

The algorithm suggested by Yan and Fan in [32] corresponds in the above taxonomy to the 
pseudo-Newton algorithm. We see that for obtaining a least-squares solution to the nonlinear 
inverse problem y = g(z), if g is holomorphic, then the Yan and Fan suggestion can result in a 
good approximation to the Newton algorithm. However, for nonholomorphic least-squares inverse 
problems and for other types of optimization problems (including the problem considered by Yan 
and Fan in [32]), the approximation suggested by Yan and Fan is not guaranteed to provide a good 
approximation to the Newton algorithm. However, as we have discussed, it does result in an 
admissible generalized gradient descent method in its own right, and, as such, one can judge the 
resulting algorithm on its own merits and in comparison with other competitor algorithms. 



Equality Constraints. The classical approach to incorporating equality constraints into the prob- 
lem of optimizing a scalar cost function is via the method of Lagrange multipliers. The theory of 
Lagrange multipliers is well-posed when the objective function and constraints are real-valued 
functions of real unknown variables. Note that a vector of p complex equality constraint condi- 
tions, 

g(z) = G 

^"^Such a claim might he true. However, it would have to be justified. 
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is equivalent to 2p real equality constraints corresponding to the conditions 

Reg(z) = G and Img(z) = G . 

Thus, given the problem of optimizing a real scalar-valued loss function i(z) subject to a vector 
of p complex equality constraints h(z) = 0, one can construct a well-defined lagrangian as 

£ = £(z) + Re g(z) + Xj Im g(z) , (147) 

for real-valued j>-dimensional lagrange multiplier vectors Xr and Xj. 
If we define the complex lagrange multiplier vector X by 

A = A^ + jA, GCf 

it is straightforward to show that the lagrangian (147) can be equivalently written as 

i: = £(z)+ReA-^g(z). (148) 

One can now apply the multivariate CM-Calculus developed in this note to find a stationary 
solution to the Lagrangian (148). Of course, subtle issues involving the application of the z, c- 
complex, and c-real perspectives to the problem will likely arise on a case-by-case basis. 

Final Comments on the 2nd Order Analysis. It is evident that the analysis of second-order 
properties of a real- valued function on C" is much more complicated than in the purely real case, 
perhaps even dauntingly so. Thus, it is perhaps not surprising that very little analysis of second 
properties can be found in any single location in the literature. By far, the most illuminating is 
the paper by van den Bos [25], which, unfortunately, is very sparse in its explanation.^^ A careful 
reading of van den Bos indicates that he is fully aware that there are two interpretations of c, viz the 
real interpretation and the complex interpretation. This is a key insight. As we have seen above, 
it provides a very powerful analysis and algorithm development tool which allows us to switch 
between the c-real interpretation (which enables us to use the tools and insights of real analysis) 
and the c-complex perspective (which is shorthand for working at the algorithm implementation 
level of z and z). The now-classic paper by Brandwood [14] presents a development of the complex 
vector calculus using the c-complex perspective which, although adequate for the development of 
first-order algorithms, presents greater difficulties when used as a tool for second order algorithm 
development. In this note, we've exploited the insights provided by van den Bos [25] to perform 
a more careful analysis of second-order Newton and Gauss-Newton algorithms. Of course, much 
work remains to explore the analytical, structural, numerical, and implementation properties of 
these, and other second order, algorithms. 



That I could find. Please alert me to any relevant survey references that I am ignorant of. 
'Likely a result of page limitations imposed by the publisher. 
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7 Applications 



1. A Simple "Nonlinear" Least Squares Problem - 1. This is a simple, but interesting, problem 
which is nonlinear in z G C yet linear in c G C C C^. 

Let 2 G C be an unknown scalar complex quantity we wish to estimate from multiple iid noisy 
measurements, 

yk = s + nk, 

k = 1, • • • , n, of a scalar signal s G C which is related to z via 

s = 9iz), g{z) = az + /3z. 

where a G C and /3 G C are known complex numbers. It is assumed that the measurement noise Uk 
is iid and (complex) Gaussian, Uk ~ iV(0, cr^/), with cr^ known. Note that the function g{z) is both 
nonlinear in z (because complex conjugation is a nonlinear operation on z) and nonholomorphic 
(nonanalytic in z). However, because the problem must be linear in the underlying real space 
7^ = (a fact which shows up in the obvious fact that the function g is linear in c), we expect 
that this problem should be exactly solvable, as will be shown to indeed be the case. 

Under the above assumptions the maximum likelihood estimate (MLE) is found by minimizing 
the loss function [15]^' 



1 " 

fc=i 

1 

= -y^AlVk- 0!Z - Pz\\^ 
k=l 

I " 

" l^'^iVk- - pz){yk- az - I3z) 



k=l 



^ n 

^ ^(yk -oiz- (3z){yk -az- (3z). 



k=l 

Note that this is a nonlinear least-squares problem as the function (^(z) is nonlinear in z!^^ Further- 
more, g{z) is nonholomorphic (nonanalytic in z). Note, however, that although g{z) is nonlinear 
in z, it is linear in c = {z, z)'^, and that as a consequence the loss function i(z) = i(c) has an exact 
second order expansion in c of the form (92), which can be verified by a simple expansion of i{z) 
in terms of z and z (see below). The corresponding c-complex Hessian matrix (to be computed 
below) does not have zero off-diagonal entries, which shows that a loss function being quadratic 
does not alone ensure that Hzz = 0, a fact which contradicts the claim made in [32]. 



^^The additional overall factor of - has been added for convenience. 

n 

Recall that complex conjugation is a nonlinear operation. 
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Defining tlie sample average of n samples {,^1, ■ ■ ■ ,C,k}^y 



Ik 

n — ' 

the loss function £(2) can be expanded and rewritten as 

2l{z) = (^\y\^) + aPz^ -{a{y)+B{y))z+{\a\^ + \P\^)z-z-{a{y) + p{y))z + aPz^ (149) 

or 

Since this expansion is done using the z-perspective, we expect that it corresponds to a second 
order expansion about the value z = 0, 

i{z) = m + + ^c^7^L(0)c (150) 

with 
and 

''''^ ' 2\ 2a(3 \af + \f3[ 
And indeed this turns out to be the case. Simple differentiation of (149) yields, 

^^ = aPz+^{\af + \pf)z-^ {a {y) + (3 {y)) 

^ = ap-z+^{\af + \pf)z-^ {a {y) + P {y)) 

which evaluated at zero give the linear term in the quadratic loss function, and further differentia- 
tions yield, 

which is independent of z. Note that, as expected. 



dl{z) _ d£{z) 
dz dz 

If we set the two partial derivatives to zero, we obtain two stationarity equations for the two 
stationary quantities z and z. Solving for z then yields the least-squares estimate of 2;,^^ 

^op' = I |2 ^ 1^,2 (" (y) - (y)) ■ 

\a\ - \p\ 



99 



Note that this answer reduces to the obvious solutions for the two special cases a = and /3 = 0. 
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This solution can also be obtained by completing the square on (150) to obtain 

Copt (^cc 



An obvious necessary condition for the least-squares solution to exist is that 

\af^\P\\ 

The solution will be a global'"" minimum if the Hessian matrix is positive definite. This will be 
true if the two leading principal minors are strictly positive, which is true if and only if, again, 
\af 7^ \(3f. Thus, if \af ^ the solution given above is a global minimum to the least squares 
problem. 

The condition = \pf corresponds to loss of identifiability of the model 

g{z) = az + [3z . 

To see this, first note that to identify a complex number is equivalent to identifying both the real 
and imaginary parts of the number. If either of them is unidentifiable, then so is the number. 

Now note that the condition \a\^ = \(3f says that a and j3 have the same magnitude, but, in 
general, a different phase. If we call the phase difference 0, then the condition \af = \f3f is 
equivalent to the condition 

which yields 



g{z) = e^'f'pz + (3z = e^^(3 {e^^z + e'^^z^ = e^^ (3 [e^^ z + e^2z^ = e^^f3Re je^'^^} . 



Thus, it is evident that the imaginary part of e^^ z is unidentifiable, and thus the complex number 
e^^ z itself is unidentifiable. And, since 

z = e~^^ (^e^^z^ = e~^^ ^Reje-' ^^;! + j Im je-' 2 2;| j , 

it is obvious that z is unidentifiable. 

Note for the simplest case of a = /?(</) = 0), we have 

g{z) = az + az = aRe {z} 

in which case Im {z}, and hence z, is unidentifiable. 



100 



Because the Hessian is independent of z. 
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2. A Simple "Nonlinear" Least Squares Problem - II. The "nonlinearity" encountered in the 
previous example, is in a sense "bogus" and is not a nonlinearity at all, at least when viewed from 
the c-real perspective. Not surprisingly then, we were able to compute an exact solution. Here, 
we will briefly look at the Newton and Gauss-Newton algorithms applied to the simple problem of 
Example 1. 

In the previous example, we computed the Newton Hessian of the least-squares loss function 
(149). The difference between the Newton and Gauss-Newton algorithm resides in the difference 
between the Newton Hessian and the Gauss-Newton Hessian. To compute the Gauss-Newton 
Hessian, note that 

2/ = (7(c) = (a/3)Q =Gc 
and therefore (since the problem is linear in c) we have the not surprising result that 

GAc = ^Ac 
dc 

with 

G={a(3). 

In this example, the least-squares weighting matrix isW = I and we have 

which is seen to be independent of c. From (122), we construct the Gauss-Newton Hessian as 



n^^ =p{G G) = = -(^ |„|2^|^|2)=w, 



cc 



showing that for this simple example the Newton and Gauss-Newton Hessians are the same, and 
therefore the Newton and Gauss-Newton algorithms are identical. As seen from Equations (130) 
and (132), this is a consequence of the fact that g{c) is linear in c as then the matrix of second 
partial derivatives of g required to compute the difference between the Newton and Gauss-Newton 
algorithms vanishes 

d_l'dg_^ 
dc \dc^ 

From the derivatives computed in the previous example, we can compute j as 



■^-f*)** all 



'dfm" ({"^)"\ 




dc ) \lm£l\" ] UsM")"; 2V 2a/3 |a|- + |/J|V V*/ 



K. Kreutz-Delgado — Copyright © 2003-2009, All Rights Reserved - Version UCSD-ECE275CG-S2009vl.O 70 



or 

/C 



C = C„pj - c . 



The optimal update in the Newton algorithm is therefore given by 

Ac.-,K,r.(M^)^-,«„-.(^)'' 

The update step in the Newton algorithm is given by 

c„e„ = c + aAc . 
If we take the "Newton stepsize" a = 1, we obtain 

c„,„ = c + Ac = c + c„p, - c = c„p, 

showing that we can attain the optimal solution in only one update step. For the real case, it 
is well-known that the Newton algorithm attains the optimum in one step for a quadratic loss 
function. Thus our result is not surprising given that the problem is a linear least-squares problem 
in c. 

Note that the off-diagonal elements of the constant-valued Hessian H^^ are never zero and 
generally are not small relative to the size of the diagonal elements of TY^c- This contradicts the 
statement made in [32] that for a quadratic loss function, the diagonal elements must be zero.^°^ 
However, the pseudo-Newton algorithm proposed in [32] will converge to the correct solution when 
applied to our problem, but at a slower convergent rate than the full Newton algorithm, which is 
seen to be capable of providing one-step convergence. We have a trade off between complexity 
(the less complex pseudo-Newton algorithm versus the more complex Newton algorithm) versus 
speed of convergence (the slower converging pseudo-Newton algorithm versus the fast Newton 
algorithm). 



3. The Complex LMS Algorithm. Consider the problem of determining the complex vector 
parameter a G C" which minimizes the following generalization of the loss function (2) to the 
vector parameter case, 

£(a) = E {|efc|^} , ek=r]k-a^^k, (151) 

for ?7fc G C and .^^ G C". We will assume throughout that the parameter space is Euclidean so that 
= I. The cogradient of £{a) with respect to the unknown parameter vector a is given by 

""it is true, as we noted above, that for the quadratic loss function associated with a holomorphic nonhnear inverse 
problem the off-diagonal elements of the Hessian are zero. However, the statement is not true in general. 
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To determine the cogradient of 

\ekf = ekCk = CkCk = {r]k - a^^k){Vk - a^^k] 



note that 



Cfc = iVk - a^ik) = iVk - iko) 
and that = (r/fc — cb^ik) is independent of a. Then we have 

= -ekg-Jka 



The gradient of \ek\ = CkCk is given by 



/ (9 \ ^ 



Thus, we readily have that the gradient (direction of steepest ascent) of the loss function i{a) = 
E{|efc|^} is 

V, iia) = -E {^kCk} = -E {^k iVk - a) } . 

If we set this (or the cogradient) equal to zero to determine a stationary point of the loss function 
we obtain the standard Wiener- Hopf equations for the MMSE estimate of a.'"^ 

Alternatively, if we make the instantaneous stochastic-gradient approximation, 

Vai{a) ^ VaOak) = Va|efc|2 = -a-cfc = a {Vk - ^k^k) , 

where is a current estimate of the MMSE value of a and — Va^(a) gives the direction of steepest 
descent of i{a), we obtain the standard LMS on-line stochastic gradient-descent algorithm for 
learning an estimate of the complex vector a, 

ttk+i = ak- akVJ^ak) 

= ak + akikek 

= ak + akik {Vk - ik'^k) 

= {I - ak^k^k) ak + ak^kVk ■ 

Thus, we have easily derived the complex LMS algorithm. 

Complex LMS Algorithm: ak+i = (/ - c^k^k^k) «fc + «fc6% • (152) 



'"-Which, as mentioned earher, can also be obtained from the orthogonaHty principle or completing the square. 
Thus, if the Wiener-Hopf equations are our only goal there is no need to discuss complex derivatives at all. It is only 
when a direction of steepest descent is needed in order to implement an on-line adaptive descent-like algorithm that 
the need for the extended or conjugate derivative arises. 
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